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I.  INTRODUCTION 

This  document  is  s final  contract  report  for  contract  DAAD05- 
74-C-0749.  Tha  research  under  this  contract  Is  aimed  toward  the 
development  of  a computer  code  which  utilizes  the  method  of  charac- 
teristics to  solve  the  problem  of  two-phase  flow  with  shocks.  This 
type  of  flow  Is  typical  of  flow  in  a projectile  launch  tube,  where 
combustion  products  and  unburnt  propellant  grains  are  mixed. 

Usually,  this  complicated  flow,  which  may  contain  shocks  of  relatively 
large  amplitude,  Is  approximated  either  by  a pure  gas  flow  or  by  a 
simplified  two-phase  flow  model  where  the  gas  and  particles  have  the 
same  or  predetermined  relative  velocities.  A better  model  to  handle 
this  problem  Is  required  to  test  the  accuracy  of  simplified  formu- 
lations and  to  produce  a more  accurate  prediction  of  the  flow, 
particularly  in  the  presence  of  shocks. 

In  addition  to  presenting  the  general  formulation,  a dis- 
cussion of  the  proper  initial-boundary  values  to  be  prescribed  for 
the  present  mixed  hyperbolic-parabolic  equations  is  presented.  It 
is  shown  that  there  is  a region  immediately  behind  the  bullet  in 
which  only  the  gas  phase  is  present.  This  consideration  requires 
the  development  of  an  interface  between  a one-phase  and  two-phase 
region.  The  computer  code  1WDFL0  which  utilizes  the  method  of  charac- 
teristics and  finite  difference  techniques  to  solve  the  problem  of  a 
projectile  accelerating  in  a gun  barrel  is  presented. 
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1. 


State  of  the  Art 


Quite  a few  textbooks  and  nanuals,  such  as.  Corner  [1],  Hunt  [2] 
and  the  Amy  Design  Handbook  [3]  exist  In  the  field  of  internal 
ballistics;  however,  most  of  these  works  give  only  qualitative  de- 
scriptions, or  present  simple  experimental  results.  They  are  not 
sufficient  for  modern  design  application;  more  un-to-date  computer 
codes  are  needed  to  better  modal  the  internal  ballistics  problem. 

In  recent  years,  computer  codes  have  been  developed  to  treat 
the  Internal  ballistics  problem.  Examples  of  which  are  codes  written 
by  Baer  and  Frankie  [4]  and  Baer  [5].  These  codes  for  the  most  part 
are  limited  in  that  they  are  restricted  to  one-phase  flows  where  the 
particle  motion  Is  predetermined.  This  is  satisfactory  when  the 
particles  are  very  small,  and  the  loading  ratio  is  low,  thus  the 
particles  and  the  gas  are  in  equilibrium.  However,  for  problems  where 
finite  size  propellant  particles  are  packed  with  a high  loading  ratio, 
the  motion  of  the  particles  is  different  from  that  of  the  gas  and  the 
flow  should  be  treated  as  two  distinct  phases.  In  1956,  under  BRL 
sponsorship,  a group  at  the  University  of  Maryland  discussed  the 
equations  of  two-phase  flow  and  proposed  schemes  for  numerical  calcu- 
lation; however  no  attempt  was  made  to  write  a code  [6].  Other  works 
such  as  Refs.  [7]  and  [8]  are  only  for  special  applications,  with 
various  limitations  and  simplifications. 

There  is  a need  for  an  up-to-date  two-phase  flow  code  that  will 
trace  shock  waves  exactly,  handle  high  loading  ratios,  include  the 
effects  of  finite  particle  volume,  accommodate  a general  form  of  the 
equation  of  state  and  various  other  features.  With  today's  high  speed 
computers  and  the  knowledge  of  numerical  methods,  it  is  feasible  to 
develop  such  a code. 

To  solve  this  problem  numerically  either  a two-dimensional  or  a 
one-dimensional  code  may  be  developed.  As  discussed  by  Moretti  [9] 
two-dimensional  codes  are  certainly  more  versatile;  however,  for 
certain  problems,  the  results  of  two-dimensional  codes  are  often  less 
than  satisfactory,  with  uncertainties  in  convergence,  stability  and 
physical  meaning.  On  the  other  hand,  one- dimensional  codes  can  be 
used  more  conveniently  for  studying  the  Importance  of  different  para- 
meters and  for  delineating  the  physical  nature  of  the  problem.  In  view 
of  the  present  state-of-the-art  in  two-phase  flow,  a one  dimensional  code 
would  be  most  useful  to  solve  interior  ballistic  problems. 

In  this  report,  we  present  a one-dimensional  two-phase  flow 
code,  in  which  a shock  wave  is  traced  exactly  both  in  front  of  and 
behind  the  projectile.  In  continuous  regions,  we  shall  use  the  basic 
method  of  characteristics  supplemented  by  finite-difference  techniques 
in  places  where  characteristics  do  not  exist.  The  method  of  character- 
istics has  the  advantage  that  it  reveals  more  directly  the  flow  prop- 
erties and  wave  structure. 
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2.  General  Features  ot  tht  Coda 


The  following  are  some  ot  the  feature!  that  are  included  in  the 
final  code. 

a.  The  computer  code  ia  written  in  aucn  a way  that  any  consistent 
set  of  units  can  be  uaea. 

b.  ihe  code  i?<  written  in  FORTRAN  IV  language. 

c.  Provisions  nave  seen  made  to  account  for  neat  loss  through  the 
barrel,  wall  friction,  and  masa  loaa  through  holes  in  the 
barrel. 

d.  The  code  is  designed  to  handle  gradual  changes  in  bore  area. 

e.  The  code  does  not  handle  ignition;  therefore,  all  propellant 
particles  will  be  assumed  to  be  Ignited  before  calculations 
begin.  An  ignition  code  must  be  used  to  generate  the  initial 
conditions  p(x,0),  e(x,0),  p(x,0),  pp(x,0),  z(x,0),  u(x,0) 

and  u^(x,0)  which  will  be  substituted  into  the  present  code. 

A shock  wave  may  be  present  initially. 

f.  The  bore  resistance  is  treated  as  a function  of  x,  u and  Ap. 

g.  The  gas  in  front  of  the  bullet  is  treated  including  the 
possible  formation  of  a shock. 

h.  A single  shock  wave  can  automatically  be  inserted  and  traced 
behind  the  projectile. 

i.  The  barrel  configuration  is  of  the  form  shown  in  Fig.  (1). 


Figure  1 - Schematic  Barrel  Configuration 
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II.  GENERAL  FORMULATION 

1.  Lif  ratura  Review 

The  problem  of  two-phase  flow  has  drawn  considerable  attention  In 
recent  yeare.  This  Interest  la  stimulated  mostly  by  applications 
such  as  solid  propellant  rockets,  nuclear  raactors,  fuel  sprays, 
lunar  ash  flow  and,  of  course,  interior  ballistics.  The  governing 
equations  of  two-phase  flow  have  been  presented  by  many  authors.  A 
few  of  these  will  be  discussed  hare. 

Rudinger  [10,  Section  4]  derived  a set  of  equations  for  solid  parti- 
cles suspended  in  a gas  where  the  particle  volume  fraction  is  finite. 

A simpler  set  of  equations  neglecting  the  particle  volume  was  also 
presented  [10,  Section  7].  This  simple  set  of  equations  was  solved 
numerically  for  an  unsteady  flow  problem  [11].  Rudinger  gave  detailed 
physical  meaning  to  various  terms  in  his  equations  and  discussed  the 
shock  conditions  and  relaxation  times.  His  formulation  Is  not  suitable 
for  direct  application  to  the  present  problem,  because  (a)  it  is  limited 
to  an  ideal  gas  equation  of  stats,  (b)  no  mass  transfer  between  the 
particle  and  the  gas  is  considered,  (c)  no  transport  of  mass,  momentum, 
and  energy  between  the  mixture  and  the  duct  wall  is  allowed,  (d)  no 
provision  is  made  for  a change  of  duct  area  along  the  axis,  and  (e) 
the  pressure  gradient  force  on  the  particle  is  neglected. 

Migdal  and  Agosta  [12]  derived  two-phase  flow  equations  by  consider- 
ing the  particle  terms  as  a source  of  drag  and  heat  transfer  in  a 
pure  gas  flow.  They  included  the  effects  of  mass,  momentum  and  energy 
transport  between  the  particles  and  the  gas,  without  giving  explicit 
expressions  for  them  and  limited  their  study  to  the  case  of  small 
particle  volume.  Using  concepts  in  continuum  mechanics,  such  as  par- 
tial stress  and  partial  energy,  Soo  [13]  derived  equations  for  multi- 
phase flow.  Hr  emphasized  the  importance  of  the  size  distribution  of 
the  solid  particles,  and  wrote  a set  of  fundamental  equations  treat- 
ing solid  particles  of  different  sizes  as  distinct  species  in  the 
mixture.  This  approach  is  not  directly  suitable  for  practical  appli- 
cation because  in  general  the  size  distribution  of  the  solid  particles 
is  more  or  leaa  continuous.  Panton  [14]  treated  the  two— phase  flow 
problem  by  defining  all  physical  quantities  in  terms  of  their  time 
and  special  averages.  This  approach  may  be  suitable  for  a turbulent 
flow  study,  but  is  too  cumbersome  for  other  purposes.  Murray  [15]  de- 
rived a set  of  equations  to  be  uved  primarily  for  fluidization  appli- 
cations. He  assumed  a constant  particle  volume  ratio  and  neglected 
pressure  forces  acting  on  the  particles.  In  each  of  his  equations  of 
motion  terms  Involving  the  time  derivative  of  both  the  fluid  and 
particle  velocities  are  included.  As  will  be  shewn  later,  these  terms 
will  change  the  nature  of  the  governing  equations.  Marble  [16]  also 
applied  modern  techniques  of  fluid  mechanics  to  the  two-phase  flow 
problem  of  a gas  and  solid  particles,  however,  he  limited  his  study 
to  the  caee  of  a negligibly  small  particle  volume  fraction. 
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Most  studies  of  a two-phase  flow  are  limited  to  the  case  of  a 
negligibly  small  particle  volume  fraction.  As  mentioned  before, 

Rudlnger  [10]  Included  the  volume  fraction  terms,  but  did  not  consider 
the  pressure  gradient  force  acting  on  the  particle.  Pai  [17]  Included 
the  particle  volume  fraction  term  and  also  the  pressure  gradient  force. 

He  treated  this  pressure  force  by  considering  the  particles  as  a pseudo- 
fluid  with  a partial  pressure.  This  pseudo-fluid,  however,  does  not 
contribute  to  the  pressure  of  the  mixture.  In  doing  this,  the  term 
containing  the  spaclal  derivative  of  the  particle  volume  fraction  , — , 

d X 

was  not  included.  For  problems  where  e is  not  negligible  the  de- 
rivative may  be  of  the  same  order  as  other  terms  in  the  equations  and 
its  omission  may  cause  appreciable  error. 

In  this  report,  we  derive  the  governing  equations  including  the 
pressure  gradient  force  on  the  particles,  and  also  the  term.  It  will 

3x 

be  shown  that  the  inclusion  of  this  term  causes  the  characteristic  prop- 
erties of  the  equations  to  change.  When  this  term  is  neglected,  two 
compatibility  equations  exist  along  the  triple  degenerate  characteristic, 
" up,  while  when  it  is  included,  only  one  exists  (for  further  details 

see  Sec.  111,1) 

2.  Basic  Assumptions 

The  governing  equations  are  derived  on  the  basis  of  the  following 
assumptions: 

a.  The  equation  of  state  of  the  gas  is  of  the  form  p-p  (p,E)  or 
p-p  (p,T). 

b.  The  average  size  of  the  actual  particles  will  be  used  in  the 
equations;  they  are  incompressible;  their  specific  heat  is 
constant,  and  the  temperature  is  uniform  within  each  particle. 

c.  The  particles  are  uniformly  distributed  over  the  cross-section 
of  the  duct,  and  their  size  and  average  spacing  are  small  com- 
pared with  the  cross-sectional  area. 

d.  The  flow  is  treated  as  one- dimensional,  thus  changes  in  the 
cross-sectional  area  of  the  duct  must  be  sufficiently  gradual. 

Abrupt  changes  in  the  cross-sectional  area  must  be  treated  by 
matching  the  continuous  flow  regions  on  both  sides. 

e.  The  drag  force  between  the  gas  and  particle  phases  is  assumed  known. 
It  may  be  prescribed  as  a function  of  any  of  the  flow  variables. 

The  modified  Stokes  formula  and  Ingebo  formula,  which  were  used 

by  Pal  and  Rudlnger,  are  typical  drag  force  expressions. 

f.  The  effect  of  the  particles  on  the  gas  flow  is  distributed  over 
the  entire  gas  phase  by  mixing.  This  mixing  involves  only  a 
small  gas  volume  and  la  therefore  assumed  to  take  place 

ins  tantaneous ly . 
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g.  The  list  of  the  particles  to  be  considered  will  be  a few  orders 
of  magnitude  larger  than  the  molecules  of  the  gas.  It  will  be 
assumed,  therefore,  that  the  particles  do  not  contribute  to  the 
pressure  of  the  mixture.  The  pressure  of  the  gas-pertlcle 
mixture  Is  given  by  the  pressure  In  the  gas  phase  alone.  The 
volume  fraction  of  the  particles,  e,  will  not  be  assumed  smell. 
Therefore,  the  pressure  gradient  will  act  on  the  particles,  as 
well  as  on  the  gas,  and  will  be  Included  In  the  momentum  equa- 
tions of  the  particles. 

h.  Gravitation  and  other  body  forces  will  not  be  Included  in  the 
equations. 

1.  The  density  ratio  between  the  particles  and  the  gas  is  small 
enough  so  that  terms  containing  p/p  may  be  dropped  from  the 
equations.  p 

j.  The  particle  density  pp  Is  constant. 

3.  Governing  Equations 

Under  these  assumptions,  the  governing  equations  for  the 
gas  medium  are  (for  a detailed  derivation  see  Appendix  A)  the  continuity 
equation. 


Do  . 3u  cm  dA  , 

— + a t — ■ — — -r-  + to  - u» 
Dt  3x  A dx  w 


the  momentum  equation, 


5?  + ^ H - - ?t<F-v  - »<V">  + “«(vu» 


o 3x  o 
and  the  energy  equation 

|E  £|u_.eu|£<i1  [Q4<j  + U(F-F  ) - Fu  - u(uu  - j u2  + E) 
Dt  p 3x  o 3x  o v w'  p P 2 


(1) 


(2) 


1 2 

+ a)  (uu.  - -z  u + E)  - u>  „ . 

w w 2 wp  Pp  A 3x 


_E_  _ Pu^r.cJL  M I 

« A * 


(3) 


The  governing  equations  for  the  particles  are  the  continuity  equation 


DPo 


3u 


a u .. 

. _ P P M _ 
A dx 


u - w 


wp 


Dt  + °p  3x 
the  momentum  equation, 

V * {FF  1(F-Fw>  * “<VU)  + 

+ <F  + V ‘ WV’ 

and  the  energy  equation,  which  is  uncoupled  from  the  system, 

DP1  _ 3u  pu  . . 

— -2.  + -E — _E  + — E — - — [q  + o 
Dt  p^  3x  3x  o„  iyp  ^wp 

P P P 

.2 


(4) 


(5) 


u-  ep  ur 

+ u F_  + (E  + -f  ) (urH*) ) - 


p WP 


wp' 


M j 
dx  1 


(6) 
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where  D/Dt  and  D?/Dt  are  material  derivatives  of  the  gas  and 
particle  properties,  respectively,  and  are  given  by 


D 

_±  + 

3 

Dt  " 

at  + 

U 3x 

DP  . 

9 

3 

Dt  " 

+ 

°P  3x 

Equations  (l)-(  6)  represent  a system  of  6 equations  in  terms  of 
8 unknowns  a,  u,  E,  p,  o^,  u^,  E^  and  e.  To  complete  this  system 

we  must  supply  the  equation  of  state  of  the  gas 

P - P (P,E)  (.  7 ) 

and  the  definitions  of  o and  a 

P 

o - (l-e)p  ( 8 ) 

and 

o ■ c P„  ( 9 ) 

P P 

With  the  assumption  that  Pp  is  constant  (the  particle  is  incom- 
pressible), Eqs.  ( 1) — ( 9 ) represent  a system  of  9 equations  in 
terms  of  9 unknowns  o,  u,  E,  p,  p,  e,  o„.  Up,  and  Ep.  The  terms 
“*  V wwp»  F*  Fw*  V Q*  V V and  Vp  ®PPearin8  on  the  right 
hand  side  of  Eqs.  (1)  to  (9)  will  be  considered  as  given  input  infor- 
mation; they  may  be  either  functions  of  x and  t,  or  functions  of  the 
flow  properties.  Usually,  they  do  not  contain  derivatives  of  the  flow 
properties  with  respect  to  x or  t;  therefore  they  will  not  affect  the 
characteristic  directions  and  the  form  of  compatibility  equations 
associated  with  the  system. 


t 

! 


4.  Initial  and  Boundary  Condi t Iona 

When  attempting  to  determine  a set  of  properly  posed  initial  and 
boundary  conditions  for  a set  of  equations,  it  is  beneficial  to  have 
a thorough  understanding  of  their  characteristic  directions  and  com- 
patibility equations,  since  the  two  are  closely  related.  The  details 
concerning  the  characteristics  associated  with  Eqs.  (1) — ( 5)  can  be 
found  in  Sec. (Ill, 1);  however,  & few  pertinent  observations  will  be 
presented  here. 

First,  the  particle  phase  equations,  Eqs.  (4  ) and  ( 5).  are  weakly 
coupled  to  the  gas  phase  equations,  Eqs.  (1),  (2)  and  (3)  (note  that 
the  converse  is  not  true),  thus  the  characteristic  directions  associated 
with  each  phase  can  be  calculated  separately.  Second,  the  compatibility 
equations  associated  with  Eqs.  (4)  and  ( 5 ) can  be  calculated  inde- 
pendent of  Eqs.  (1),  (2)  and  (3),  however,  again,  the  converse  is  not 
true.  The  full  set  of  equations,  Eqs.  (l)-( 5 ) must  be  used  to  calcu- 
late the  compatibility  equations  for  the  gas  phase. 

With  this  information  and  the  understanding  that  boundary  con- 
ditions are  strongly  dependent  on  the  characteristic  equations,  we  shall 
make  the  assumption,  that  in  determining  the  proper  boundary  conditions 
for  our  problem,  we  may  treat  the  three  gas  equations  and  the  two 
particle  equations  separately.  The  gas  equations  are  then  treated  as 
completely  hyperbolic,  and  possess  the  same  boundary  conditions  as 
one-phase  compressible  flow.  The  one  phase  equations  incorporating  an 
ideal  gas  assumption  will  now  be  presented  briefly.  They  are 

3t+“^+p  al*° 
la  + u in  + I it  . o 

3t  3x  p 3x 

JL  (JL)+U  — (-2- 

3t  V PY  3x  ' PY  ) - 0 


where 


n r 2y/(y-1> 

£-•  (f-) 

p0  c0 

The  characteristic  directions  and  their  corresponding  compatibility 
equations  for  this  system  are  then 


i dx 
along  - u, 


i dx 
along  ^ 


u + c. 


and  along  ^ ■ u - c, 


d(  ) - 0 
PY 

du  + -^-dp  ■ 0 

du ^-dp  » 0 

YP  r 


For  hyperbolic  systems  a properly  posed  set  of  initial  and  boundary 
conditions  is  relatively  easy  to  determine.  Courant  and  Friedrichs  [18] 
discussed  this  problem  in  detail  while  introducing  the  concept  of 
"space-like"  and  "time-like"  curves.  Essentially,  they  simply  state 


r 


that  for  each  characteristic  reaching  a point  on  a boundary  from  outside 
the  domain,  one  dependent  variable  must  be  specified.  As  an  example,  let 
us  take  a typical  boundary  where  the  boundary  coincides  with  a particle 
path,  dx/dt  ■ u.  If  this  is  a right  boundary,  the  characteristic  dx/dt  - u-c 
falls  outside  the  region  of  interest,  while  if  it  is  a left  boundary  the 
characteristic  dx/dt  * u+c  falls  outside.  In  either  case  one  dependent 
variable  must  be  specified  on  each  boundary. 

In  determining  the  proper  initial  conditions,  we  treat  the  initial  time 
line  as  a boundary  defined  by  t*0.  For  this  case  all  three  characteristics 
fall  outside  the  domain  and  thus  all  three  dependent  variables,  p,  p and  u 
must  be  specified. 

For  the  problem  presented  in  this  report,  we  specify  three  gas  proper- 
ties on  the  initial  time  line  and  one  gas  property  on  each  boundary  as  can 
be  seen  in  the  summary  at  the  end  of  this  section. 

The  determination  of  proper  boundary  conditions  for  the  particle 
equations  can  not  easily  be  placed  on  a rigorous  foundation  since  they  are 
parabolic  in  nature.  It  is  felt,  therefore,  that  insight  into  this  problem 
could  be  gained  by  studying  a few  classical  parabolic  systems  where  proper 
boundary  conditions  have  been  established.  Based  on  observation  of  these 
cases,  we  shall  propose  a hypothesis  on  boundary  conditions,  and  apply  it 
to  our  parabolic  particle  equations. 


Boundary  Condition  Hypothesis 

For  two  first  order  partial  differential  equations  in  terms  of  two 
dependent  variables  u and  v,  and  two  independent  variables  £ and  n, 

1.  If  dn  ■ 0 is  a degenerate  characteristic,  then  there  can  be  only  one 
boundary  at  n * constant,  along  which  either  u or  v may  be  prescribed. 

No  boundary  conditions  need  be  prescribed  on  any  other  n * constant  line. 

2.  If  one  of  the  independent  variables  is  time  t,  then  only  one  boundary 
with  t - constant  can. exist,  and  further, 

a.  If  t - constant  is  a degenerated  characteristic,  either 
u or  v may  be  specified; 

b.  If  t ■ constant  is  not  a characteristic,  both  u and  v 
must  be  specified,  which  is  a typical  Cauchy  initial 
value  problem. 


Consider  the  system  of  equations 
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which  Is  equivalent  to  one  second  order  equation 


• a • 


where  the  omitted  terms  do  not  contain  any  derivatives.  For  our  present 
purpoie,  we  shall  use  the  first  order  system,  Instead  of  the  single  second 
order  equation.  It  can  be  shown  that  for  (10),  dn  * 0 is  a characteristic 
line.  According  to  our  hypothesis,  the  boundary  of  the  domain  must  be 
"open"  in  the  positive  n direction.  Boundary  condition  can  be  specified 
along  only  one  n - constant  line  ([19],  p.  692).  If  £ is  not  the  time 
coordinate,  then  two  boundaries  exist  along  two  £ ■ constant  lines 

A typical  problem  of  this  nature  is  the  heat  transfer  problem  governed 
by  the  equation 

32T  2 3T 

3x2  " * 9t 

or  the  equivalent  first  order  system 

3G  2 3T 
3x  " ® 3t 


where  T is  temperature  and  a is  a constant.  The  characteristic  direction 
associated  with  this  equation  is 

dt  - 0 . 

A properly  posed  set  of  initial  and  boundary  conditions  for  this  problem, 
as  shown  in  Fig.  2,  is 

t - 0 specify  T 

x - 0 specify  T 

x ■ L specify  T 

where  L is  the  thickness  of  the  plate. 


T Specified  on  All  Three  Boundaries 
Figure  2 Boundary  Conditions  for  the  Hast  Transfer  Problem 
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Another  example  is  the  boundary  layer  equations.  If  only  the  terms 
relevent  to  our  present  discussion  are  Included,  the  governing  equation 
may  be  written  as 


Sy2 


1 3u 

V 


which  has  an  equivalent  first  order  system 


3G  1 3u  , 

__  . — u — + • • • 

3y  v 3x 


where  x Is  in  the  direction  of  flow  and  lies  in  the  plane  of  the  plate,  y 
is  normal  to  the  plane  of  the  plate,  u is  the  velocity  in  the  x direction 
and  v Is  the  kinematic  viscosity.  The  characteristic  direction  for  this 
system  is 

dx  - 0 

A properly  posed  set  of  boundary  conditions  for  this  problem  as 
shown  in  Fig.  3 is 

y ■ 0 specify  u 

y “>  specify  u 

x - xQ  specify  u (xQ,y) 


Figure  3 Boundary  Conditions  for  Boundary  Layer  Flow 


In  both  of  these  examples,  £ is  not  the  time  coordinate.  In  our 
particle  phase  case,  the  equations  may  be  represented  by 
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DPo  3u 


DPu 


Dt 


where  the  independent  variable  equivalent  of  £ of  Eq.  (10)  la  time. 
Therefore,  only  one  £ - constant  boundary  exists,  as  shown  in  Fig.  4. 


u or  o 
P P 
specified 


^-Boundary  does  not  exist 
J (2nd  t ■ constant  line) 


, ✓'Boundary  does  not  exist 
tj  (x  ■ constant  ia  a 
characteristic) 


W * w * 


> r * * -r-r-rr-T  r r r ? 

o and  u specified 
P P 

(t  - constant  is  not  a characteristic) 


♦ x 


Figure  4 Boundary  Conditions  for  the  Particle  Phase  of  Our  Flow 


In  the  actual  problem,  the  right  hand  side  of  the  domain  is  bounded  by 
a gas-particle  interface,  see  Fig.  S.  No  boundary  conditions  need  be 
specified  on  this  boundary  for  the  particles,  the  solution  of  the  gas  equations 
and  the  particle  equations  will  yield  the  location  of  this  interface  line. 


Figure  5 Physical  Plane  Description  of  TWo  Phase  Flow  Behind  a Projectile 
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In  sutraary,  the  boundary  conditions  for  our  problem  are  Breech  (x*0) 
u(0,t)  ■ 0 
Up(0,t)  - 0 


Left  Bullet  surface  (x  ■ Xg)  and  (x  < Xj^) 


“b  dt 
duB  1 

dt"  " " ^FR_(pBEF“PAFT)AB^ 
u(xB,t)  - Ug 

Right  Bullet  surface  (x  - Xg  + Lfi) 

u(Xg  + Lg,t)  - u(Xg,t) 

Muzzle  - after  bullet  has  left  (x  ■ Xj^) 

P ^XM,t>  " PBEX  + ^PM~PBEX^  ^^EX^BD 

when  [u(xM,t)  < c(xM,t)] 
or 

u(xM,t)  - c(xM,t) 

Muzzle  - before  bullet  has  left  (x  ■ Xj^) 

p(xM,t)  - pM  when  [u^.t)  < c(xM,t)] 
or 

uCXjj.t)  - c^.t) 
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Muzzle  - after  shock  has  left  (x-x^) 

pC^.t)  - pmT  + (PM"PExu^t“tEXlT^tSD 

when  [u(xM,t)  < cCs^t)! 
or 

u(xM,t)  - c(^j,t) 


Initial  time  line  (t**0) 

u(x#0)  - f^(x) 

p(*»0)  - f2(x) 

P<x,0)  - f3(x) 

z<*iO)  - f^(x) 

Up(x,0)  - f5(x) 

Pp(x,0)  • constant 
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5. 


Shock  Waves 


The  treatment  of  shock  waves  presented  here  is  based  on  the 
hypothesis  that  when  a shock  wave  passes  through  a particle-laden  gas, 
the  changes  occur  at  such  a rapid  rate  that  immediately  behind  the 
wave  the  particles  have  not  had  sufficient  time  to  react.  The  shock 
wave  produces  a rapid  deceleration  of  the  gas  accompanied  by  a nearly 
discontinuous  rise  in  pressure.  Therefore,  as  presented  by  Krlebel  [21] 
and  Rudinger  [10],  the  initial  disturbance  caused  by  the  shock  wave 
is  not  influenced  by  the  presence  of  particles  in  the  gas.  All 
properties  behind  the  shock  can  be  calculated  from  the  shock  conditions 
if  the  gas  properties  in  front  of  the  wave  are  known  along  with  the 
shock  speed  or  one  gas  property  behind  the  wave.  This  state  immed- 
iately behind  the  wave  is  known  as  the  "frozen"  state. 


The  equations  used  to  calculate  the  gas  properties 
state  are  the  standard  Rankine-Hugoniot  shock  relations 

in  the  frozen 

(U-up  - p2  (U-u2) 

(11) 

P2_P1  “ P1  (U"ui>  (U2"’V 

(12) 

(13) 

where  subscripts  1 and  2 represent  properties  in  front  of  and  behind 
the  shock,  respectively,  and  U is  the  shock  speed.  As  previously 
mentioned,  the  properties  of  the  particles  are  identical  in  front  of 
and  behind  the  wave. 

6.  Equation  of  State 

In  formulating  the  governing  equations,  it  was  only  natural  to 
assume  that  the  equation  of  state  was  a function  of  the  specific 
Internal  energy,  E,  and  the  density,  p or 

P ■ p(p,E)  (14) 

However,  in  the  practical  application  of  flow  in  a gun  barrel,  the 
equation  of  state  will  more  than  likely  be  given  in  the  form 

p » p'  (p,T)  (15) 

therefore,  provision  must  be  made  to  Incorporate  this  form  into  our 
formulation.  In  this  section  we  will  treat  the  subject  in  detail 
and  present  specific  examples. 

From  the  mathematical  standpoint,  if  we  can  arrive  at  a re- 
lation between  E,  T,  and  p of  the  form 

f(E,T,p)  - 0 (16) 

then  Eqs.  (IS)  and  (16)  will  combine  to  be  equivalent  to  Eq.  (14). 


The  procedure  used  to  generate  Eq.  (16)  involves  the  solution  of 
several  partial  differential  equations  from  classical  thermodynamics 
[22].  In  order  to  solve  these  equations  two  conditions  must  be  speci- 
fied. First,  Cv(pq,T)  must  be  given,  which  for  this  problem  it  is 
assumed  to  be  of  the  form 

Cv(p0'T)  " A1  + A2  T 


Secondly,  the  constant  E(Pq,Tq)  must  be  given.  Knowing  these  two 
conditions,  we  can  proceed  to  determine  Eq.  (16). 


We  begin  by  calculating  the  constant  volume  specific  heat  of 
the  gas  by  solving  the  equation 


Then,  we  calculate  E(p,T)  and  subsequently  Eq.  (16)  from  the  relations 


3E 

3T 


C 


v 


and 


3E 

dp 


P 


and  the  condition  E(Pq,Tq) . 

Determining  Eq.  (16)  only  conceptually  solves  the  problem  of 
handling  an  equation  of  state  in  the  form  of  Eq.  (15).  At  this  point 
we  will  treat  the  actual  technique  for  implementing  Eqs.  (is)  and  (i&) 
into  the  computer  code. 

When  the  value  of  p must  be  calculated,  assuming  that  p and  E 
are  known,  Eq.  (16)  is  solved  numerically  for  T and  then  T and  p are 
substituted  into  Eq.  (15)  to  yield  p.  When  the  derivatives  |P  and  |P. 

are  required,  they  are  calculated  from  Eqa.  (15)  and  (k^  through  the 
following  formulas 

IE.  ISL'  i£fej£L 

3E  3T  3E 

i£  . i£'  3TCP.J.EJ  + l£' 

3p  3T  3p  3p 

3T 

Where  the  value  of  -r=-  and  r—  are  found  from  Eq.  (16)  using  the  relations 

at,  dp 

il  . if  , 3f 

3E  " ’ 3E  ' 3T 

and 

3T  . . if  , 3f 

3p  3p  ' 3T 
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We  will  now  present  two  equations  of  state  which  have  been 
implemented  Into  the  computer  code.  The  first  is  the  "virial  equation" 
given  by 

p - P ^ T [1  + p B(T)  + p2  C(T) ) (17) 

and  the  second  is  the  van  der  Waala  equation  given  by 

<■  ■ ’ s 1 < ° °2  (18> 

where  R/M,  n,  and  a are  constants  which  are  material  dependent  and 
B(T)  and  C(T)  are  empirically  determined  functions  of  temperature. 

It  should  be  noted  that  if  a is  set  equal  to  zero  Eq.  (18)  is  the 
Noble-Abel  equation  of  state  and  if  both  a and  n are  set  equal  to 
zero  it  becomes  the  ideal  gas  equation  of  state. 

The  values  of  Cv  and  f(p,E,T)  corresponding  to  Eq.  (17)  are 
given  by 

2 

Cy(p,T)  - \ + A2T  - (p-pQ)  |L(I  [B(T) 

3T 

+ j (p+pQ)  C(T) ] } 


f(p,E,T) 


- E - E(pq.T0)  - A1(T-Tq)  - j A2(T2-T2) 

[ + 1 <e+co)  1 ■ 0 


while  the  values  corresponding  to  Eq.  (18)  are  given  by 
cv  <».T)  - (»0.T)  ■ A1  + A2T 

f(o.E.T)  > E - E(Tq(0q)  - Ajd-Ig)  - | Aj  (T2-!2) 

+ o(p-Pq)  - o 

As  additional  equations  of  state  become  available  they  may  easily  be 
incorporated  into  the  code. 
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7.  Area  Discontinuity  (gas-particle  interface) 

In  order  to  examine  the  behavior  of  a two-phase  material  as  it 
passes  through  an  area  discontinuity,  moving  at  a velocity  up,  let  us 
construct  a control  volume  as  shown  in  Fig.  6. (it  should  be  noted  that  for 
the  gas  particle  interface  up  ■ il).  It  is  convenient  to  treat  section  B as 
an  exit  from  the  "main"  flow  field  and  section  A as  an  entrance  back 


Section  B 


Figure  6.  Area  Discontinuity 

into  the  field.  Thus  section  A may  be  considered  as  the  entrance  of 
region  II,  and  section  B as  the  exit  from  region  I.  These  two  regions 
are  coupled  by  constraints  relating  their  respective  properties  across 
the  control  volume,  region  III.  In  the  limit,  we  may  let  region  III 
shrink  to  a line  located  at  the  area  discontinuity.  We  will  now 
determine  these  constraints.  First  we  note  in  Fig.  7 that  four  of  the 
characteristic  lines  fall  outside  the  main  flow  field:  namely  char- 
acteristics 4B,  at  section  B,  and  1A,  2A  and  3A  at  section  A.  Equations 
written  along  these  lines  must  therefore  be  replaced.  In  Sec.  (111,1) 
it  is  shown  that  one  compatibility  equation  is  associated  with  each  of 
the  four  above-mentioned  characteristic  directions;  in  addition,  the 
particle  continuity  equation  is  written  in  finite-difference  form  along 
the  characteristic  3A  (d|  _ y^).  Thus  there  are  five  missing  equations 

which  must  be  replaced  by  five  new  relations  coupling  the  properties 
of  section  A with  those  of  section  B. 

Let  us  now  determine  these  new  equations.  Since  the  particle 
equations  are  weakly  coupled  to  the  gas  equations,  the  characteristic 
directions  can  be  associated  distinctly  with  either  the  gas  or  the 
particle  phase.  We  shall  assume  that  the  equations  written  along  "gas" 
dx 

characteristics  (~  - u + c,  u)  should  be  replaced  by  equations  relating 

gas  properties,  and  that  equations  written  along  the  "particle"  char- 
acteristic (&  • Up)  should  be  replaced  by  equations  relating  particle 
properties.  dt  Therefore  we  will  need  three  equations  governing  the 
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Figure  7.  Physical  plane  (x,t)  plot  of  area 

discontinuity  showing  characteristic 
network. 
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change  In  gas  properties  between  sections  A and  B,  and  two  equations 
relating  particle  properties. 

First,  let  us  consider  the  relations  for  the  gas  phase.  Our 
treatment  essentially  follows  that  given  by  Shapiro  [23].  In  his 
approach,  mass  and  energy  are  conserved  between  the  two  sections; 
momentum  conservation  is  not  considered,  because  to  do  so  would  re- 
quire knowledge  of  the  forces  acting  between  the  wall  and  the 
fluid.  In  place  of  momentum  conservation,  Shapiro  makes  the  simpli- 
fying assumption  that  entropy  is  continuous  across  the  area  change. 
For  our  problem,  it  seems  that  this  assumption  is  satisfactory.  How- 
ever, if  entropy  calculations  are  prohibitively  difficult  because  of 
the  complexity  of  the  particular  equation  of  state  chosen,  then 
either  the  assumption  of  constant  temperature  between  the  sections 
or  the  momentum  equation  utilizing  an  approximation  for  the  pressure 
on  the  wall  may  be  used. 


Let  us  now  determine  the  two  equations  governing  the  change  in 
particle  properties^  One  relation  which  Is  easily  obtained  is  con- 
servation of  mass  between  the  sections,  namely: 

V(Waa  " Vv  ' “d)ab 

The  remaining  equation  is  more  difficult  to  determine.  As  was  the 
case  for  the  gas  phase,  applying  the  principle  of  conservation  of 
momentum  would  require  knowledge  of  the  Interaction  force  between 
the  wall  and  the  particles;  thus  rendering  its  use  impractical.  The 
equation  of  energy  conservation  for  the  particle  phase,  although  it 
could  be  applied,  would  not  complete  the  system  of  equations  since 
it  is  uncoupled  from  the  system,  and  thus  would  merely  add  one  equa- 
tion and  one  new  variable  Ep.  Entropy  calculations  are  impossible, 
for  the  model  ignores  the  thermodynamic  aspects  of  the  particle  phase. 


One  possible  approach  is  to  assume  that  the  particle  velocity 
does  not  vary  across  the  area  discontinuity.  This  supposition  is 
reasonable  when  one  considers  that,  although  the  gas  phase  velocity 
increases  instantaneously,  the  particles,  being  more  massive,  do  not 
accelerate  as  quickly.  (Such  an  approximation  brings  to  mind  the 
assumption  of  frozen  flow  used  to  calculate  property  changes  across 
shock  waves;  see  section  11(5).)  For  area  discontinuities,  this 
approximation  seems  valid  if  t is  small  or  if  the  area  change  is  from 
small  to  large;  however,  inconsistencies  (such  as  e>l)  may  arise  with 
flows  in  which  a highly  concentrated  particle  phase  travels  from  a 
large  area  to  a smaller  one.  For  this  reason,  we  have  chosen  to  write 
the  last  of  the  five  equations  in  the  form: 


(u  ).  “ K (u  )Rt 
p A m p B 

where  K is  a factor  which,  most  likely,  depends  on  geometry  and  flow 
m 

parameters  and  must  be  determined  experimentally. 


* Note  that  for  the  gas-particle  Interface  these  equations  are  not 
necessary  since  particles  do  not  exist  at  Section  B. 
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For  the  present,  we  have  made  the  convenient  assumption  that 


which  when  combined  with 


A B’ 

the  mass  conservation  equation  yields 

K ■ -r— 
m Aa 


It  is  yet  to  be  seen  how  the  factor  K affects  the  overall  flow  cal- 
culations. To  summarize,  the  equations  which  will  be  used  are,  for 
the  gas  phase. 


Conservation  of  mass 

°AiuK-%)kA  " °B(VUD)AB 

Conservation  of  energy 


« + t>a-  ch  + t>, 


Continuity  of  entropy 


SA-SB 


and  for  the  particle  phase, 
Conservation  of  mass 

0A(V“d)AA  • °B^Up,-UD)AB 
Particle  velocity  assumption 


<V* 


W» 


(19) 

(20) 
(21) 

(22) 

(23) 


8.  Maas  Transfer  Between  Solid  and  Gas  Phase 


Following  the  procedure  used  to  calculate  the 
released  while  the  propellant  is  burning  we  define 
individual  particle,  M , as 


M - 

P 

The  rate  at  which  the  mass  of 
burning,  ft  , is  then  given  by 


ft  - 
P 


V a 
P P 

a single  particle  is 


Dt  p 


D*v 


■ P 


p Dt 


rate  at  which  heat  is 
the  mass  of  an 


changing  due  to 


where 

DPVp  x DPz 

Dt  ’ S (z)  Dt"~ 
P 

and  Sp(z)  is  the  burning  surface 


21 


i 

% 


S 

\ 

! 


i 


Defining  Che  number  of  particles  per 

N . -L 

v p v 

p p p 


unit  volume  of  mixture,  N,  as 


we  then  can  define  the  rate  at  which  mass  is  being  added  to  the  gas 
phase  per  unit  volume  of  mixture,  u,  as 


<u 


- ft  N 
P 


(24) 


9.  Burning  Law  (Regression  Equation) 


Before  introducing  the  specific  equation  specifying  the  regression 
rate,  let  us  first  define  the  regression  distance  Z(x,t).  To  do  this 
we  define  a parameter  to  be  a characteristic  dimension  of.  a propellant 
grain  in  the  sense  that  it  is  the  least  dimension  which  has  to  be 
traveled  by  the  burning  surface  in  order  to  bum  the  propellant  cost* 
pletely.  The  regression  distance  is  then  defined  as  the  amount  that  p 
has  decreased  from  its  initial  value  at  a particular  x and  t.  The 
regression  rate  is  then  simply  DPZ/DT.  The  specific  form  currently  being 
used  is  the  non-linear  burning  law  given  by 


DpZ 

Dt 


a(T)  p* 


(25) 


10.  Heat  Released  During  the  Burning  of  the  Propellant 

In  order  to  calculate  q,  the  amount  of  heat  released  during  a given 
time  Interval,  we  must  be  able  to  calculate  the  volume  of  the  propellant 
burnt  at  any  instant,  V (x,t);  or  equivalently,  the  instantaneous  particle 

volume  V (x.t)  recognizing  that  the  two  are  related  we  can  write 
P 

DPV*'t>_  (26) 

Dt  " “ Dt 

If  we  now  define  a parameter  K_  as  the  amount  of  heat  released  per  unit 
volume  of  propellant  burnt,  then  thi  rate  at  which  an  individual  particle 
la  giving  off  heat  is  thus  given  by 


(27) 


Now  defining  N(,  the  number  of  particles  per  unit  mass  of  gas  as 


Mg  d o V_  p 


P P 


a V 


(28) 
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we  can  calculate,  q,  the  rate  at  which  heat  is 
phase  per  unit  mass  of  gas  by 


q - 


being  added  to  the  gas 


(29) 


Equation  (29) 


can  be  combined  with  Eqs.  (24), 
DPV 

« ■ - ",  kq  tit  ■ rV 

P 


(27)  and  (28)  to  yield 


(30) 


11.  Calculation  of  the  Initiation  and  Tracing  of  a Shock  Wave 
Behind  the  Bullet 

Here  we  present  techniques  used  to  insert  and  trace  a shock  wave 
behind  the  bullet.  These  proceedures  will  handle  shock  reflection 
from  both  the  breach  and  the  bullet  precisely;  however,  some  order  of 
approximation  will  be  necessary  to  treat  a shock  passing  through  a 
gas-particle  Interface.  It  is  proposed  to  cal- 

culate the  singularity  occurring  at  this  interaction  point  precisely 
and  then  "smear"  the  reflected  wave  to  simplify  future  calculations. 

The  wave  structure  of  this  interaction  point  is  shown  in  Figs.  (8) 
and  (9).  Figure  (8)  shows  the  interaction  of  a right- traveling 
shock  with  the  gas-particle  interface.  The  proceedure  for  calcu- 
lating this  singularity  is: 

a.  Calculate  the  properties  in  regions  a,b  and  c and  assign  them  to 
the  mesh  points  as  shown. 

b.  Calculate  the  properties  in  regions  2,3  and  5 (region  4 notation 
is  not  used  due  to  programming  considerations)  (a<->l  and  c<->6) 

c.  Assign  the  properties  in  regions  2,3,5  and  6 to  mesh  points  as 
shown.  This  has  the  effect  of  causing  the  properties,  starting 
at  point  I,  to  vary  gradually  until  they  reach  the  values  of 
region  2 just  before  the  gas  particle  interface.  A precise 
treatment  would  cause  the  properties  in  region  1 to  be  gradually 
reached  at  the  gas  particle  interface  and  then  change  instan- 
taneously to  those  in  region  2,  see  Fig.  (10).  (The  properties  in 
region  1 are  not  retained  for  future  calculation). 

It  should  also  be  noted  that  the  contact  surface  is  not  traced 
in  subsequent  calculations. 

Figure  (9)  shows  the  interaction  of  a left  traveling  shock  with  the 
gas-particle  Interface.  The  calculation  procedure  here  is  quite 
similar  to  that  used  for  a right  traveling  wave  and  is  as  follows: 
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Figure  8.  A Right  Traveling  Shock  Wave  Interacting 
With  a Gas-Particle  Interface. 
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Figure  9.  A Left  Traveling  Shock  Wave  Interacting 
With  a Gas-Particle  Interface. 
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Figure  10.  Pressure  vs  Position  Plot  Showing  Singularity 
Properties  Before  and  After  Simplification 
for  a Right  Traveling  Shock  Intersecting  with 
a Gas-Particle  Interface. 


Figure  11.  Pressure  vs  Position  Plot  Showing  Singularity 
Properties  Before  and  After  Simplification 
for  a Left  Traveling  Shock  Intersecting  with 
a Gas-Particle  Interface. 


a.  Calculate  the  properties  in  regions  a,b  and  r and  assign  them 
to  the  mesh  points  as  shown. 

b.  Calculate  the  properties  in  regions  2,3  and  5 (region  4 notation 

is  not  used  due  to  programming  considerations),  (a  1 and  c 6) 

c.  Assign  the  properties  in  regions  1,2,3  and  6 to  mesh  points  as 
shown.  This  has  the  effect  of  causing  the  properties,  starting 
at  point  II,  to  vary  gradually  until  they  reach  the  values  of 
region  3 just  after  the  Interface.  A precise  treatment  would 
cause  the  properties  in  region  6 to  be  reached  just  after  the 
Interface  and  then  change  Instantaneously  in  two  steps  to  those 
in  regions  5 and  3,  see  Pig.  (11).  (The  properties  in  regions  5 
and  6 are  not  retained  for  future  calculations) 

It  should  be  noted  that  the  contact  surface  is  not  traced  in 
subsequent  calculations. 

The  simplifications  that  we  use  in  treating  the  shock-interface 
singularity  and  in  the  subsequent  calculations  seem  to  produce  errors 
which  are,  for  the  most  part,  dependent  on  the  mesh  size  (the  smaller 
the  mesh,  the  more  accurate  the  smoothing)  and  the  magnitude  of  c 
(the  smaller  the  value  of  e,  the  smaller  the  reflected  wave).  One 
will  have  to  judge  on  an  individual-problem  basis  the  magnitude  of  these 
errors. 


III.  NUMERICAL  PROCEDURE  (TOE  METHOD  OF  CHARACTERISTICS) 

1.  Calculation  of  the  Characteristic  Directions  and 
Compatibility  Equations 

Before  calculating  the  compatibility  equations  corresponding  to 
Eqs.  (1)  to  (11)  let  us  eliminate  the  partial  derivatives  of  p and  c 
from  Eqe.  (2),  (3)  and  ( 6 ) by  using  Eqs.  ( 7),  ( 8)  and  (9  ).  We 
then  arrive  at  the  partial  differential  equations  that  are  used  in 
the  computer  code,  namely: 

gas  continuity 


O.  + U 0,  + 

ou,  ■ GC 

(31) 

t X 

X 

gas  momentum 

u,.  + u u,  + 

Ac,  +BE,  + D 0 -GM 

(32) 

t X 

X X p,x 

gas  energy 

E, . + u E,  + 

F u,  + G 0 - GE 

(33) 

t X 

X p,x 

particle  continuity 

0 . + u 0 

+ 0 u - PC 

(34) 

p,t  p p,x 

P P.X 

particle  momentum 

u . + u u 

- PM 

(35) 

p,t  p p,x 

particle  energy 

E _ + u E 

+ H u + J a «PE 

(36) 

p,t  p p,x 

P.X  P,x 

. 1 

A ■ — p,  ; 
0 r,p  * 

(1-e) 

*'  9 P>E  ' 

D " pp  u-«> p’ 

; r • f ; 

P P 

P 0 

P P 


GC,  GM,  GE,  PC,  PM,  and  PE  are  the  right  hand  side  of  Eqs.  (1)  to  (11), 
respectively  and  the  notation  X,y  represents  the  partial  derivative  of 
X with  respect  to  y. 
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The  characteristic  directions  and  compatibility  equations  for 
this  system  of  equations  has  been  calculated  using  both  the  directional 
derivative  approach  and  the  determinant  approach  producing  identical 
results.  Only  the  determinant  approach  will  be  presented  here. 

In  applying  the  determinant  approach,  we  treat  the  time  and  spa- 
tial derivatives  of  the  dependent  variables  as  unknown  quantities. 

If  we  add  the  six  continuity  equations 


do 

■ 

o,x  dx  + o,t  dt 

(37) 

du 

- 

u,x  dx  + u,t  dt 

(38) 

dE 

- 

E,  dx  + E.  dt 
’x  t 

(39) 

do 

m 

o dx  + a dt 

(40) 

P 

P.x  Ptt 

du 

m 

u dx  + u dt 

(41) 

P 

P.x  P.t 

dE 

m 

E dx  + E dt 

(42) 

P 

P.x  p,t 

to  our  system  of  partial  differential  equations  we  arrive  at  a 
system  of  12  equations  in  terms  of  12  derivatives  which  when  written 
in  matrix  notation  becomes: 


Let  N be  the  determinant  of  the  coefficient  matrix  appearing  In 
Eq.  (43)  and  be  the  determinant  of  the  matrix  formed  from  the 
coefficient  matrix  with  the  1th  column  replaced  by  the  column  vector 
on  the  right  hand  aide  of  Eq.  (43).  The  aolutlon  for  the  derlvatlvee 
may  then  be  written  In  the  form 

M*o,x  - Mj  ; N*o,t  - M2  ; etc*  (44) 

The  characterlatlc  direction*  for  thle  ayatem  of  equation*  are 
defined  aa  direction*  in  the  x,  t-plane  which  cauae 


In  aolving  for  these^directlona,  It  la  convenient  to  reduce 
the  12  x 12  determinant  N through  column  operation*  and  Laplace 
expansion  [24],  to  the  following  form, 


N 


dx-udt 

*o  dt 

0 

0 

0 

0 

-A  dt 

dx-udt 

-B  dt 

-D  dt 

0 

0 

0 

-F  dt 

dx-udt 

-G  dt 

0 

0 

0 

0 

0 

dx-udt 

-c i dt 
P 

0 

0 

0 

0 

0 

dx-u  dt 
P 

0 

-J  dt 


-H  dt  dx-udt 


(45) 


dx-udt 

-o  dt 

0 

dx-iyit 

-o  dt 

0 

-A  dt 

dx-udt 

-B  dt 

X 

0 

dx-Updt 

0 

0 

-F  dt 

dx-udt 

-J  dt 

-H  dt 

dx-iyit 

N*«  (dx  - udt)[dx  - (u  ♦ c)dt][dx  - (u  - c)dt]  (dx  - vyit)*  ■ 0 

where 

«*  -oA  +(.)(«  -£  + 

P 
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Thus  the  characteristic  directions  are 


and 


dx 

dt 

dx 

dt 


dx 

Tt  ‘ u+  c 


dx 

dt 


u - c 


u 


(46) 


dx 


where  — ■ up  is  a triple  degenerate  root. 

It  can  be  seen  that  in  order  to  have  a finite  solution  for  the 
partial  derivatives  in  Gq.  (44),  the  must  equal  zero  when  N 
equals  zero.  The  setting  of  the  determinants  equal  to  zero  leads 
either  to  an  identity,  zero  equals  zero,  or  to  one  compatibility 
equation  corresponding  to  each  of  the  characteristic  directions 
(note  that  it  is  possible  for  as  many  as  three  distinct  compatibility 
equations  to  exist  along  the  triple  degenerate  root  - Up).  For 
this  system,  only  one  distinct  compatibility  equatiofrexists  along 
•jjjr  * Up,  which  indicates  that  the  system  is  not  totally  hyperbolic 
and  thus  other  equations  in  addition  to  the  compatibility  equations 
are  needed  to  complete  the  solution. 

The  compatibility  equations  for  this  system  are 

i dx 
along  ~r—  ■ u 

2 F 

G[o  du  - (u  -u)do  ] + (u  -u)  (dE do) 

P P P P P o 


- {G[op  PM  - (up-u)  PC]  + (up-u)2  [GE  - £ GC]}dt 
dx 


(47) 


along  ■ u + c 


(G  x B + D x c)  [o  du  - (u  -u  + c)  do  ] 

- P P P P 

+ (up-u  + c)2  [B  dE  + cdu  + Ado]* 

'(G  x B + D x c)  [o  PM  - (u  -u  + c)  PC] 

— p p 

+ (up-u  + c)2  [B  x GE  + c x GM  + A GC]}dt 
dx 


(48) 

t'^9) 


and  along  • up 


du 


PM  dt 


(50) 


As  mentioned  before,  there  are  two  compatibility  equations 
missing.  Thus,  two  additional  equations  must  be  supplied.  We  will 
follow  a procedure  similar  to  that  used  in  Refs.  [6]  and  [9] 
for  one-dimensional  two-phase  flow  and  Refs.  [25]  and  [26]  for 
two-dimensional  flows  and  write  Eqs.  (34)  and  (3&)  in  finite  dif- 
ference form  along  the  particle  path  line  of  the  particle  phase 
(dx/dt  * Up)  namely:  31 


PC 


(51) 


and 


o 

P 


p,x 


Hu  + J o 
P.x  p,x 


PE 


(52) 


The  addition  of  Eqs.  (51)  and  (52)  to  our  numerical  procedure 
necessitates  the  calculation  of  Up|X  and  0p,x  ftt  the  new  point. 

To  accomplish  this  we  follow  a procedure  similar  to  that  used  in 

Ref.  [6]  and  write  the  "continuity"  equation  for  Up  and  op  along 

the  gas  characteristics,  dx/dt  - u + c,  namely: 

• 

d u ■ [ U (u  + c)  + u ] dt  (53) 

P P.x  p,t 

and 

do-  [ o (u  + c)  + o ] dt  (54) 

P P.x  — p,t 

Equations  (47)  to  (54),  after  being  written  in  a second  order 
accurate  finite-difference  form,  will  be  used  to  generate  a 
solution  to  the  problem. 

Our  general  numerical  procedure  will  utilize  the  Hartree 
(constant  time)  scheme.  In  this  scheme,  it  is  assumed  that  the 
values  of  all  dependent  variables  are  known  at  discrete  mesh  points 
lying  on  a constant  time  line.  A new  time  line  is  established  to 
meet  a stability  criterion  and  the  properties  are  calculated  at 
points  where  the  gas  particle  path  originating  from  the  known  points 
on  the  old  time  line  intersects  the  new  time  line. 

The  only  problem  tl.  t arises  in  the  use  of  this  technique 
concerns  Eqs.  (53)  and  i54).  As  can  be  seen,  the  time  and  spacial 
derivatives  of  Up  and  ap  must  be  available  on  the  old  time  line  be- 
fore the  new  time  plane  can  be  calculated.  If  the  values  of  Up  and 
op  and  not  their  derivatives  are  specified  on  an  initial  constant 
time  line  the  x derivatives  will  be  calculated  using  3 point  forward 
or  backward  difference  schemes  for  the  left  and  right  boundaries 
respectively  and  central  difference  for  Interior  points  and  then 
Eqs.  ( 4.)  and  (5)  will  be  solved  for  the  time  derivatives.  These 
values  will  then  be  stored  for  use  in  calculating  the  next  time  line. 


2,  Mathematical  Implications  of  the  Physical  Assumptions 

Before  examining  the  effects  that  our  physical  assumptions 
have  on  the  characteristics,  let  us  see  what  we  can  conclude  in 
general  about  our  system  of  equations.  To  do  this  let  us  partition 
the  determinant  N?  Eq.  (45),  into  4 minors,  namely 
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where 


Oi 

! o2 

N * 

... 

1 . . . 
| 

Oj 

i 04 
1 

dx-udt 

-a  dt 

0 

0i  ■ 

-A  dt 

dx-udt 

-B  dt 

0 

-F  dt 

dx-udt 

0 

0 

0 

02  - 

-D  dt 

0 

0 

-G  dt 

0 

0 

0 

0 

0 

03  - 

0 

0 

0 

0 

0 

0 

l 

i 


I 


^4 


dx-Updt 

0 

-J  dt 


-o  dt 
dx-Updt 


-H  dt 


0 

0 

dx-Updt 


Applying  Laplace  expansion  to  N we  can  see  that  if  all  the  terns 
in  Q3  are  zero  (the  particle  equations  are  weakly  coupled  to  the 
system,  [24]). 

N*-  IqJ  X |Qj  . 

Thus,  as  long  as  Q3  - 0,  any  change  to  the  terms  in  Q2  will  not 
effect  the  characteristic  directions  of  the  system. 
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Now  let  us  examine  the  effects  of  a)  ignoring  the  terms  con- 
taining e,x  and  b)  Including  terms  which  contain  a factor  p/Pp. 

Case  (a)  is  equivalent  to  setting  D,  G and  J equal  to  zero. 
Immediately  one  can  see  that  setting  D and  G equal  to  zero  cannot 
effect  the  characteristic  directions  of  the  system  because  they 
appear  only  in  Q2.  Examination  of  Q4  quickly  shows  that  J does 
not  enter  into  the  evaluation  of  the  determinant  and  therefore  Its 
value  doesn't  effect  the  characteristic  directions.  The  effect  of 
neglecting  e,x  on  the  compatibility  equations  is  not  quite  as 
simple  to  see;  however,  after  expanding  the  M*  determinants  one 
finds  that  setting  this  term  equal  to  zero  yields  two  distinct 
compatibility  equations  along  dx  ■ Up. 

The  Inclusion  of  terms  containing  p/Pp,  case  (b) , has*  the 
effect  of  altering  Q^,  Q2,  Q3  and  by  adding  a term  Involving  the 
acceleration  of  the  gas  to  the  particle  momentum  equation,  a term 
involving  the  acceleration  of  the  particle  to  the  gas  momentum 
equation,  and  terms  involving  both  accelerations  to  both  the  gas 
and  particle  energy  equations.  The  determinant  of  the  coefficient 
matrix,  N,  for  this  system  is  of  the  form: 


dx-udt 

-0  dt 

0 

0 

0 

0 

-A  dt 

(dx-udt) 

-B  dt 

-D  dt 

a2 (dx-Updt) 

0 

N*« 

0 

c2 (dx-udt) 

(dx-udt) 

-G  dt 

c3 (dx-Updt) 

0 

0 

0 

0 

dx-Updt 

-a  dt 

0 

! 

1 

0 

(dx-udt) 

0 

0 

b2 (dx-Updt) 

0 

0 

d2 (dx-udt) 

0 

0 

dj (dx-Updt) 

dj^  (dx-Updt) 

where 

the  coefficients  a^, 

a2’  bl»  V 

Cl»  C2’  C3 

, d.,  d.  and 

d.  are 

in  general  functions  of  the 

dependent  variables. 

As  can  be  seen. 

this  system  is  quite  a bit  more  complex  than  the  one  actually  treated 
in  this  report.  The  expansion  of  N will  result  in  a order  poly- 
nomial In  (dx/dt)^  which  must  be  solved  numerically  to  yield  the 
characteristic  directions.  The  number  of  real  characteristic  directions 
resulting  from  this  equation  is  directly  dependent  on  the  values 
of  the  coefficients;  thus,  the  nature  of  the  system,  at  lease  con- 
ceptually, could  change  as  the  values  of  the  dependent  variables 
change.  The  derivation  of  the  compatibility  equations  present  a 
similar  problem  because  they  cannot  be  determined  until  the  char- 
acteristic directions  are  known. 
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3.  Finite  Difference  Fora  of  the  Compatibility  Equations 


The  specific  finite  difference  equations  resulting  from  Eqs. 

(47)  to  (5l)  will  now  be  presented.  Referring  to  Fig.  (12)  we  see 
that  the  locations  where  the  characteristics  erainating  from  the  point 
being  calculated,  point  4,  cross  the  previous  time  line,  are  labeled 
points  A,B,C  and  D. 
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• - - - mesh  points 

0 - - - intermediate  points  (base  of  characteristics) 


Figure  12.  Characteristic  network  used  in 
finite  difference  scheme 


Utilizing  the  above  point  scheme  and  the  notation  that  XA  corresponds 
to  a property,  X,  evaluated  at  point  A,  X*A  corresponds  to  (X^)A  and 
AtA  corresponds  to  At  between  points  4 and  A:  equation  (47)  solved 
for  E4  becomes 


. 1 r_fi- 

2 L(  v"4 ) 


1 

2 


G.  0 . G„  0 „ 
4 p4  + C pC 


(“p4'u4)2  <upc'uC)2J 

1 


(up4-upc) 


+ 7 r|  (0  . - 0 ) + -r 

^pc’VJ  p4  pC  2 


F F 
-*  + -£ 
°4  °C 


(o4-°C) 


4 t. 


[ G4  gp4  PM4 

L <w2 


GC  qpC  PMC 

(upc'“c)2 


C4  PC4 
<Up4-“4> 


GC  PCC 


<“Pc-“c> 


+ GE4  + GEC  - 


F4  GC4 


FC  GCC 


(55) 
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equation  (48)  solved  for  u^  becomes 


UB‘  2 


Y.  . a . Y o 

14  p4  + IB  pB 

Y24  C4  Y2B  °B 


(u 


p4 


UPB> 


4 


Y Y 

14  + IB 


Y24C4 


Y2B  CB 


(°p4  ' °PB) 


f B. 

B 

(E.  - E_)  - | 

fA. 

A 1 

-i  + 

lc‘ 

CB 

4 B 2 

CB  J 

(°4  - °B) 


At. 


14  °p4  PM4  , Y1B  CTpB  PMB  Y14  PC4 


Y2  c 
24  4 


Y2B  cB 


Y24  C4 


Y1B  PCB  , B4  Ct4  , »B  GEB 
Y2B  CB  C4  CB 


+ G«4  + GMb  +c 


*4  0C4  . *B  GCB 


B 


(56) 


where 

and 


Yj  « G x B + D x c 
y2  ' “p  - “ ‘ c 


equatlon(49) solved  for  becomes 


°4  " °A  “ 2 


Y34  °p4  + Y3A  PpA 
Y44  A4  Y4a  AA 


W 


4 
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Y44  A4  Y4A  AA 


(o  , - o .) 
p4  pA' 


«4  - EA> 
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4 


(u4  - UA) 


At, 


Y34°P4  PM4  , Y3A°PA  PMA 


Y44  A4 


Y4A  AA 


Y34  PC4  _ Y3A  PCA  . B4  GE4 
Y44  A4  Y4A  AA  A4 


BAGEA  C4  °^4 


CA  GMA 


+ GC.  + GC. 
4 A 


where 


and 


Y^*GxB-Dxc 


Y.  * u - u + c 
4 P 


(57) 


equation  (50)  solved  for  u , becomes 

At_  P 


v * V + -r  (pm4  + PV 

(51)  solved  for  o ^ 

becomes 

At_  | 

r 

r9u  'i 

. D 

o , ■ a - + — — - 1 

pc.  + pc  - a . — E- 

p4  p0  2 

L 4 

D p4  [dx  j 

r3u  ^ 

— CT 

— E.] 

PD 

i3x  J D_ 

(58) 


(59) 


Equations  (55)  to  (59)  are  the  equations  used  to  calculate  a regular 
point  in  our  grid.  Under  certain  conditions,  however,  they  must 
be  modified. 


The  first  special  case  we  will  present  occurs  when  the  proper- 
ties at  points  A,B,C  and  D are  Identical.  When  this  occurs,  Eqs, 
(47)  to  (51)  reduce  to 

du  ■ GM  dt 
do  - GC  dt 
dE  ■ GE  dt 

d u - PM  dt 
P 

do  - PC  dt 
P 
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These  equations  are  finite  differenced  by  averaging  GM,  GC,  GE, 

PM,  and  PC  between  point  4 and  point  C. 

The  second  special  case  occurs  when  either 

u . - u.  - 0 
p4  4 

or 

u - u - 0 

pC  C 

In  this  case  the  use  of  Eq.  (55)  would  lead  to  a division  by  zero. 
Thus,  we  must  finite  difference  Eq.  (47)  In  the  following  form, 


E4  • Ec  + 


: 7 7 K V + GC  V)["T  (pm4  + ptfc> 


p4  4 pC  “C 


- <up*'"pc)]  ■ [yvV  + °c(vv] 

x (PC,  * PCC)  - ( op4  - V>]  - 


Pl«<V«-“4)' 


F1C(udC"UC) 


[-T  <GC4  + 0CC>  - <“4  - “C* 


+ -J-  (CE,  + GEC) 


4.  Stability 

He  have  conducted  an  investigation  to  gain  insight  into  the 
stability  of  a mixed  hyperbolic-parabolic  system  of  equations  such 
as  we  have  here.  Actually  deriving  a mathematically  rigorous  sta- 
bility criterion  for  the  highly  non-linear  equations  of  this  re- 
port is  quite  a formidable  task  and  not  within  the  scope  of  this 
research.  However,  we  feel  that  it  is  quite  helpful  to  have  an 
idea  of  the  type  of  problems  that  might  be  encountered  in  treating 
our  mixed  system.  It  is  our  conclusion  that  the  standard  stability 
criterion  for  hyperbolic  equations  may  not  be  sufficient  for  our 
system. 

To  demonstrate  the  above  mentioned  conclusion  let  us  look  at  a 
particular  system  of  equations  for  a problem  coupling  sound  propa- 
gation and  heat  flow  as  presented  by  Richtmeyer  [20] » namely; 


■ 


mew nef tsar  as,.  <*w'*tt* 


h - c -k '» + <*-»*> 

3w  _ 3u 
3t  " C 3x 

2 

3e  _ ^ 9 e c 3u 


where  u is  the  material  velocity,  w ■ c V/Vq,  V is  the  specific 
volume,  e ■ E/c,  E is  the  specific  Internal  energy,  c is  the 
isothermal  sound  speed,  and  o is  the  ratio  of  thermal  conductivity 
tc  specific  heat  at  constant  volume.  These  equations  are  formed 
by  coupling  the  hyperbolic  fluid  dynamic  equations  to  the  para- 
bolic heat  flow  equations.  The  first  two  equations  above  may  be 
said  to  be  hyperbolic  in  nature,  while  the  third  is  parabolic  in 
nature.  An  analogous  situation  occurs  with  the  equations  of  this 
report  where  Sqs.  (1),  (2)  and  (3)  are  hyperbolic  in  nature,  and 
Eqs.  (.4)  and  ( 5)  are  parabolic  in  nature. 

Without  going  into  the  finite  differencing  details  which  are 
presented  on  page  171  of  [20]  we  will  proceed  right  to  the  con- 
clusions concerning  the  stability.  Although  a precise  stability 

system  was  not  found.  Ref.  [20]  states 
for  it  tc  satisfy  the  stability  cri- 
fluid  dynamics  equations  and  the  un- 
The  stability  criteria  for  these  two 


< 1 


It  is  also  stated  that  in  the  limit,  as  Ax  and  At  go  to  zero,  the 
second  condition  implies  the  first  and  thus  is  assumed  to  be  the 
stability  criterion  for  the  system. 

This  example  demonstrates  the  possibility  that  a mixed  hyper- 
bolic parabolic  system  may  be  finite  differenced  in  such  a manner 
as  to  yield  a stability  criterion  which  is  directly  dependent  on 
the  parabolic  segment  of  the  equations.  With  this  in  mind  let  us 
examine  the  possibilities  that  exist  for  our  system  of  equations. 

a.  The  hyperbolic  stability  criterion  dominates  and  is  given  by 

< M + e)  K 


criterion  for  this  complete 
that  it  is  surely  necessary 
teria  of  both  the  uncoupled 
coupled  heat  flow  equation, 
systems  are  respectively: 

4 c At 


Ax 


and 


o At 
(Ax)  “ 


or 


K Ax 


t < /'I'  l-V-  \ < K, 

- ( ul  + c)  - 1 


where  K is  equal  to  1 


39 


HMMMVna 


When  this  criterion  is  programmed,  the  value  of  K is  set  somewhat 
less  than  1 to  allow  for  the  fact  that  the  properties  on  the  new 
time  line  are  unknown  when  At  is  calculated  and  must  be  estimated. 

b.  The  parabolic  stability  criterion  dominates  and  is  of  the 
same  form  as  the  hyperbolic  criterion.  This  criterion  is  programmed 
in  the  same  manner  as  the  hyperbolic  criterion  except  that  K may  now 
be  significantly  less  than  1. 

c.  The  third  possible  criterion  results  from  the  comparison 
of  the  parabolic  segment  of  our  system  with  the  heat  flow  equation 
in  Sec.  (11,4).  It  can  be  shown  that  from  a characteristic  stand- 
point the  roles  of  x and  t are  reversed  between  the  two  systems. 

The  stability  criterion  for  the  heat  flow  equation  is  of  the  form 


* 


j 

1 


*£_  < 
ix2  2 

which  leads  to  the  possibility  that  the  stability  criterion  for 
our  parabolic  equations  may  take  the  form 


At 

Ax 


> k; 


or 


At  > K2(Ax) 


1/2  . 


K. 


If  is  less  than  then  a value  of  At  may  be  chosen  such  that 
K3  < At  < K3 

However,  if  K3  is  greater  than  Ki  no  value  of  At  will  satisfy  both 
the  hyperbolic  and  parabolic  stability  criterion  and  the  system  is 
unconditionally  unstable. 

All  indications  from  past  research  are  that  this  last  case  does 
not  exist;  Rudlnger  [10]  and  Rudinger  and  Chang  [11]  have  solved 
systems  of  equations  which  are  quite  similar  to  those  presented  here 
and  have  not  reported  stability  problems.  It  is  felt  that,  although 
care  must  be  taken  in  running  the  code,  selection  of  At  to  achieve 
stability  can  be  achieved  by  using  the  stability  criterion  of  a.  and 
allowing  K to  vary  between  0 and  1. 
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5.  General  Point  Iteration  Procedure 


The  most  difficult  portion  of  the  computer  code  is  centered  on 
the  solution  of  Eqa.  (55)  to  (59).  At  first  a simple  procedure  was 
tried  where,  each  equation  is  solved  for  one  particular  variable, 
then  updated  by  averaging  the  most  recent  value  with  the  old  value 
utilizing  a relation  of  the  form 


n+1 


+ BK  x u 

aF 


n-1 


(60) 


where  the  values  of  AK  and  BK  can  be  adjusted  at  the  programmers 
discretion.  This  proved  unsuccesaful^causing  rapid  divergence  of 
the  system.  Upon  close  examination  of  the  equations,  it  was 
felt  that  the  best  hope  for  solution  would  be  to  uncouple  the  gas 
and  particle  phases  and  solve  Eqs.  (55),  (56)  and  (57)  as  a set  and 
Eqs . (58)  and  (59)  as  a set.  This  is  accomplished  by  assuming 
that  the  particle  properties  are  constant  while  solving  Eqs.  (55). 
v 56)  and  (57)  and  conversely  the  gas  properties  are  constant  when 
solving  Eqs.  (58)  and  (59).  Even  making  this  assumption,  the  solu- 
tion of  Eqs.  (55),  (56)  and  (57)  for  the  gas  properties  is  a quite 
formidable  task  due  to  the  high  degree  of  nonlinearity.  It  was 
decided  to  solve  these  three  equations  for  the  variables  , u and  E 
using  the  Newton-Raphson  technique.  Before  proceeding,  it  should  be 
noted  that  Eq.  (55)  is  linear  in  the  variable  E and  can  be  written 
in  the  form 

E - E(u,o)  (61) 

Thus  conceptually,  Eqs.  (56)  and  (57)  when  combined  with  (61)  r in 
be  written  in  the  form 


g(u,o ) - 0 
f(u,o)  - 0 


respectively.  The  Newton  Raphson  procedure  for  two  equations 
can  then  be  utilized,  namely 


and 
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u + Au 
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0 + Ao 
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Au  ■ 
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(62) 

(63) 

(64) 

(65) 
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Before  applying  this  technique  further  convenient  simplifi- 
cations were  made.  They  are  as  follows: 

a.  along  with  holding  the  particle  properties  constant  we  decided 
to  also  hold  p and  c constant.  This  assumption  was  more  or 
less  made  out  of  necessity  due  to  the  difficulty  in  evaluating 
the  partial  derivatives  of  g and  f. 

b.  Provision  was  made  and  subsequently  adopted  to  hold  the  values 
of  C.C,  GM,  GE,  PC,  PM  constant  through  this  sub  iteration. 

Having  made  these  assumptions  preliminary  calculations  produced  rapid 
convergence  of  the  Newton-Raphson  routine  and  subsequent  convergence 
of  the  entire  general  point  routine  by  then  solving  Eq.  (58y  for  Up 
and  Eq.  (59)  for  Op  and  then  correcting  the  entire  set  of  solutions 
using  equations  of  the  form  (60)  with  AK  ■ 2 and  BK  ■ 1. 


The  calculations  proceeded  routinely  for  many  time  lines  utili- 
zing this  technique;  however,  as  the  value  of  up  started  to  approach 
u (particles  became  small)  the  solution  to  the  system  of  equations 
began  to  diverge.  To  better  understand  the  cause  of  this  problem 
and  the  necessary  steps  to  correct  it  the  entire  system  of  equa- 
tions was  closely  examined.  It  was  found  that  the  trouble  was  rooted 
in  the  solution  of  Eqs.  (55),  (56)  and  (57)  and  more  specifically  in 
Eq.  (55).  To  analyze  this  problem  we  will  write  Eq.  (53)  in  a 
general  form,  namely: 


E 


f1(u,o) 

(Up-u)2 


(du  -PMdt) 
P 


f«(u,o) 

+ <^^r<dvPcdt)  + f3<u-°) 

P 


(66) 


Now  let  us  examine  a term  of  the  nature  (up-u)m  (m  ■ -1,-2)  in 
context  with  our  iteration  proceedure  noting  that  up  is  constant. 

To  do  this  it  is  simplest  to  forsake  rigor  andgo  directly  to  a 
typical  numerical  example.  First  let  us  assume  a problem  where  the 
difference  between  Up  and  u is  large;  such  as: 


u - 10 

P 

u - 100 


Now  let  us  see  what  happens  to  the  coefficient  of  the  first  term  on  the 
right  hand  side  of  Eq.  (61),  which  we  will  call  F^,  if  u varies  by, 
let  us  say  2 from  one  iteration  to  the  next,  i.e.  Au  in  Eq.  (64)  is 
equal  to  2. 


1)  for  u 


2)  for  u 


■ 100  we  have 

100 

p100  . il 

1 8100 

■ 102  we  have 


.102 


8464 
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assuming  that  f =<  f which  in  reality  is  pretty  close  to  being 

t rue  we  have  a percent  change  in  F of  approximately  4%. 

For  our  second  example  let  us  assume  a problem  where  the  dif- 
terence  between  Up  and  u is  small  such  as 


u = 100 

Now  again  let  us  assume  that  u varies  by  2 from  one  iteration  to  the 
next  and  examine  F^  or 


1)  for  u = 100 

100 

100  JL 

1 ” 25 


2)  for  u * 102 


,'.02 

1 


r102 

1 

49" 


Making  the  same  assumptions  on  as  in  the  previous  example  we  can 
see  that  the  change  in  Fj  is  now  approximately  50%.  It  turns  out 
that  this  rapid  change  in  F^  causes  the  entire  Newton-Raphson  tech- 
nique to  become  unstable. 


The  problem  now  centers  on  what  steps  should  be  taken  to  correct 
this  deficiency.  It  is  quite  obvious  that  we  cannot  tolerate  "large" 
changes  in  u during  an  Iteration  when  the  values  of  u and  up  are 
"close".  One  approach  to  the  problem  would  be  to  attempt  to  improve 
the  first  guess;  however,  first  guesses  consisting  of  the  base  point 
properties  and  the  solution  of  the  linearized  system  of  equations 
proved  unsuccessful  and  this  approach  was  abandoned. 

The  approach  that  proved  successful  is  as  follows.  First,  attempt 
to  solve  the  complete  set  of  equations.  If  this  fails,  set  fi  equal 
to  zero  and  solve  that  set  of  equations.  Then  use  its  solution, after 
having  corrected  E^  by  adding  the  term  back  in, as  a first  guess 
in  solving  the  complete  set.  If  this  still  fails,  set  both  fi  and  f2 
equal  to  zero  and  follow  the  same  proceedure. 


The  degree  to  which  the  solution  of  the  equations  with  f^  set 
equal  to  zero  approximates  the  solution  to  the  complete  set  is  re- 
lated to  how  close  the  term  (dUp  - PMdt)  is  to  zero.  Notice  that 
his  term  is  in  the  same  form  Eq.  (35)  the  only  difference  being 
that  Eq.  (61)  is  written  along  dx/dt  * u and  Eq.  (35)  along  dx/dt  - Up. 
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Thus  it  can  be  easily  seen  that  as  Up  approaches  u,  the  term 
(dtip-PMdt)  written  along  dx/dt  - u approaches  zero  and  the 
solution  of  the  simplified  system  approaches  the  solution  of  the 
complete  system.  In  summary,  it  should  be  noted  that  although  the 
term  (dup-PMdt)  is  approximately  equal  to  zero  when  the  system 
converges,  during  the  iteration  proceedure  it  can  become  quite  large 
and  cauoe  convergence  problems. 


IV.  Results  and  Discussions 


Before  presenting  a complete  test  run  for  the  M16  rifle,  we  will  briefly 
point  out  two  features  of  the  code  TWOFLO  which  are  not  treated  by  other  codes. 
The  first  feature  is  the  treatment  of  a gas  only  region  behind  the  bullet, 
typical  plot  of  e and  u vs.  time  in  the  vicinity  of  the  interface  between 
this  gas  only  region  and  the  two  phase  flow  region  is  shown  in  Fig.  (13). 


u 


Figure  13.  Plot  of  loading,  c,  and  gas  velocity,  u, 
vs.  distance  from  the  breach,  x. 


Notice  that  at  this  interface  there  Is  a discontinuity  in  c leading  to  a 
discontinuous  drop  in  particle  velocity.  The  second  feature  deals  with  a 
shock  wave  traveling  behind  the  bullet.  Figure  (14)  shows  a plot  of  the 
physical  plane  and  a plot  of  pressure  vs.  position  as  the  shock  passes 
through  the  gas-particle  interface.  Referring  to  the  pressure  curve  in 
Fig.  (14)  one  can  see  that  the  effects  of  the  contact  line  have  been  neglected 
and  the  rarefaction  wave  has  been  smoothed  out  consistent  with  the  dis- 
cussion in  Sec.  (11,7) 
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Presented  nex.  are  numerical  results  for  the  Mlb  rifle  beginning  at 
tin  time  of  complete  ignition  and  extending  until  the  bullet  leaves  the 
b.trre  1 . 

Two  separate  runs  wore  made  each  with  a different  size  propellant  particle. 
The  propellant  used  is  pancake  shaped,  (1 . 814x10“ -^Kg  of  WC  846).  The  drag 
lorce  between  the  particles  and  the  gas  is  approximated  by  that  of  a sphere 
having  the  same  radius  as  the  radius  of  the  pancake.  The  propellant  was 
assumed  to  be  compacted  by  the  initial  primer  blast.  This  is  treated  in  the 
code  by  shortening  the  rear  of  the  cartridge  by  0.005m.  The  loading  after 
compaction,  e,  was  0.575.  Other  pertinent  data  for  these  runs  can  be  found 
in  the  sample  input  section  of  Appendix  C. 

The  sizes  of  the  propellant  particle  used  for  the  runs  are  summarized 
.is  follows: 


Particle  Radius  Particle  Thickness 


RUN 

A 

2.730 

-4 

x 10  m 

3.810 

x 10 

RUN 

B 

1.365 

x 10 

1.905 

x 10 

For  Run  A the  particle  dimensions  conform  to  the  average  dimensions  of  the 
actual  propellant.  The  particle  size  of  Run  B is  one-half  of  that  of  Run  A. 
Run  B was  chosen  to  determine  the  effect  of  increasing  the  propellant  surface 
area.  The  calculated  results  for  these  two  cases  are  compared  with  the 
experimental  results  of  Trafton  (28], 

Figure  (15)  shows  the  velocity  of  the  bullet  plotted  against  the  distance 
from  the  base  of  the  cartridge.  Run  A yields  a muzzle  velocity  that  is  35% 
lower  than  the  experimental  one;  Run  B produced  a much  higher  muzzle  velocity, 
but  is  still  17%  lower  than  the  experimental  value 

Examining  Fig.  (16)  one  can  see  that  the  pressure  produced  at  the 
midpoint  of  the  chamber  by  Runs  A and  B brackets  the  experimental  values 
fairly  well,  noting  that  the  300  psec  time  duration  for  complete  ignition  to 
occur  (this  represents  the  time  between  primer  ignition  and  the  initiation 
of  TWOFLO  calculations)  is  not  exact;  any  error  in  this  time  would  shift  the 
pressure  curves  horizontally.  Although  Run  B produces  pressure  higher  than 
the  experimental  value  at  a station  in  the  chamber,  it  produces  a much  lower 
pressure  at  a station  further  downstream  in  the  barrel,  as  shown  in  Fig.  (17). 
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Figure  16.  Plot  of  pressure  vs.  time  from  ignition  at  a point  located 


Plot  of  pressure  vs.  time  from  Ignition  at  a point  ioc 
at  the  gas  port  (0.328m  from  the  base  of  the  cartridge) 


Fig.  (17)  is  a plot  of  pressure  vs.  time  at  the  gas  port, 
substantially  down  the  barrel.  By  this  time,  both  Runs  A and  B produce 
pressures  that  are  substantially  below  the  experimental  results.  These 
results  seem  to  indicate  that  too  large  of  a portion  of  the  energy,  as 
computed  by  TWOFLO,  remains  in  the  rear  section  of  the  barrel  whereas  it 
should  be  more  concentrated  near  the  bullet. 

There  are  two  possible  reasons  for  the  discrepancy  between  the 
calculated  and  experimental  results.  The  first  is  the  uncertainty  about 
the  physical  parameters  used  as  input  data  in  the  calculation.  These 
include  the  burning  rate  of  the  propellant,  the  drag  coefficient  of 
particles,  ignition  time,  the  state  at  the  end  of  the  ignition  process, 
the  friction  between  the  bullet  and  the  barrel,  etc.  Hopefully,  a 
computer  code  that  is  numerically  accurate,  together  with  certain 
easily  measured  physical  parameters,  can  be  used  to  determine  all  the 
other  parameters,  in  the  future. 

The  other  possible  reason  is  the  inaccuracy  of  the  code,  either 
in  the  governing  equations  used,  or  in  the  numerical  solution  of  these 
equations.  Preliminary  test  calculations  show  that  the  numerical 
solution  converges;  solutions  from  using  two  different  mesh  sizes 
vary  very  little.  Calculation  of  a simple  one-phase  flow  problem  also 
indicates  that  results  are  very  close  to  the  exact  solution.  Therefore, 
we  have  confidence  in  the  numerical  accuracy. 

As  discussed  in  Appendix  A,  the  governing  equations  used  involve 
certain  approximations.  In  particular,  a term  in  the  particle  momentum 
equation  has  been  neglected  for  simplicity  in  applying  the  method  of 
characteristics.  This  term  is  relatively  small  for  small  values  of  e, 
but  can  become  important  for  large  value  of  e.  It  can  be  seen  from 
Eq.  (A12)  that  by  neglecting  this  term,  we  decreased  the  particle 
acceleration.  This  might  have  contributed  to  the  slower  motion  of  the 
particles,  and  the  concentration  of  energy  and  pressure  near  the  chamber. 

For  further  development  of  this  work,  we  suggest  the  following: 

a.  Modification  of  Governing  Equations  - Although  the  equations  we  used 
are  quite  elaborate  and  include  many  more  terms  than  most  earlier 
works,  preliminary  numerical  results  indicate  that  certain 
neglected  terms  should  be  retained.  For  instance,  the  Du/Dt  term 

in  the  particle  momentum  could  be  retained. 

b.  Complete  Test  Runs  - Appropriate  sample  problems  should  be  run  with 
the  code  to  ascertain  the  convergence  and  stability  of  the  numerical 
calculation.  Comparison  with  simpler  problems  with  exact  solutions, 
and  with  other  numerical  methods  should  be  made.  A parametric 
study  to  determine  the  importance  of  various  terms  in  the 
equations,  and  the  effect  of  certain  physical  quantities  would 

also  be  desirable. 
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i.  Comparison  with  Experimental  Results  - This  type  of  comparison 
can  serve  to  determine  the  accuracy  of  the  calculated  results 
if  the  experimental  results  are  accurate,  or  vice  versa. 

t> 

r- 

d.  Solve  a Proolem  with  Shock  Waves  - This  code  lias  the  capability 

to  treat  shock  waves;  shocks  are  traced  exactly,  instead  of  being  smeared 
by  artificial  viscosity  as  in  finite-difference  methods.  The 
subroutines  are  all  debugged,  but  have  not  been  tried  out  on  a 
physical  problem  with  shocks.  This  should  be  done. 

In  conclusion,  the  following  points  about  the  TWOFLO  code  can  be 

made . 

1.  it  is  one  of  the  most  "sophisticated"  one-dimensional  codes.  It 
treats  the  two  phases  separately,  includes  the  effects  of  wall 
area  change,  wall  friction,  heat  and  mass  transfer  through  the 
wall,  includes  the  3t/Tx  effect,  etc.  The  governing  equations 
include  many  additional  terms  as  compared  with  other  existing 
codes . 

2.  It  is  accurate.  TWOFLO  incorporates  the  method  of  characteristics, 
which  is  inherently  more  accurate  than  the  finite-difference  method. 

It  handles  the  initial  and  boundary  conditions  in  a logical  manner. 


3*  It  has  the  capability  of  treating  shock  waves.  Shocks  are  traced  exactly. 
For  those  physical  problems  where  shocks  are  present,  this  code 
can  yield  more  accurate  results.  Even  in  problems  without  shocks, 
the  characteristic  lines  and  contact  lines  calculated  from  this 
code  can  reveal  more  about  the  nature  of  the  flow. 


i. 


* * 
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APPENDIX  A-DERIVATION  OF  TWO  PHASE  FLOW  GOVERNING  EQUATION’S  * 


In  deriving  the  governing  equations  for  this  system  let  us 
treat  the  gas  phase  and  the  particle  phase  separately  while  cou- 
pling their  motion  through  interaction  terms.  When  treating  each 
phase,  we  will  assume  that  it  forms  a continuum  and  occupies  the 
entire  control  volume  at  a given  instant.  In  utilizing  this  con- 
cept, we  must  replace  the  classical  density  which  represents  mass 
per  unit  volume,  with  a new  term  representing  mass  per  unit  volume 
of  mixture.  We  will  also  assume  that  the  particles  are  rigid  and 
that  they  may  be  considered  to  be  large  with  respect  to  the  molec- 
ular size  of  the  gas  and  thus  do  not  contribute  to  the  pressure 
of  the  mixture.  Having  taken  these  assumptions  into  consideration, 
we  may  write  the  governing  equation  for  each  component  of  the 
mixture. 

1.  Continuity  equations 

The  continuity  equation  for  the  gas  phase  can  be  obtained  by 
establishing  a control  volume  and  equating  the  tine  rate  of  increase 
of  mass  inside  the  control  volume  to  the  time  rate  of  mass  added 
to  the  control  volume,  where  the  mass  added  is  composed  of  three 
terms;  the  net  mass  flux  into  the  control  volume  through  the  end 
surface  normal  to  the  flow,  the  mass  addition  resulting  from  the 
particles  burning  and  the  mass  transport  through  the  wall,  or 

■ - + _ (W^A)  (Al) 


which  can  be  rearranged  as 


Dt 


3u 

3x 


au  dA 

~7~  T"  + oj  - u) 
A dx  \ 


(A2) 


Similarly,  the  continuity  equation  for  the  particle  phase  is 


at 


3(o  u A) 

P P _ 

3x 


(uiA)  - (w  A) 

WP 


(A3) 


which  can  also  be  rearranged,  producing 


DPo 


Dt 


3u  a u 

. _ p p p dA 

+ CT„  TT"  ~ . j w - <u 

P 3x  A dx  wp 


(A4) 


2.  Momentum  equations 

Before  proceeding  to  derive  the  momentum  equations  we  will 
make  several  additional  assumptions,  namely: 


We  would  like  to  make  a special  acknowledgement  to  Dr.  Alvars  Celmins  who 
made  an  invaluable  contribution  in  the  formulation  of  the  equations  in 
this  Appendix. 
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a.  When  mass  transfer  occurs  between  the  phases,  the  par- 
ticle phase  always  loses  mass,  the  gas  phase  always  gains 
mass;  or  u > 0, 

b.  The  gas  and  particle  phases  may  discharge  through  the 
wall,  but  injection  1.'  net  treated,  or 

u > 0;  u >0. 

w — wp  — 

When  leakage  occurs,  the  gas  and  particle  phases  lose 
mass  that  is  moving  ar  velocities  of  u and  u 
respectively.  One  possible  assumptionWis  to  cBnsider  the 
discharged  mass  having  the  same  axial  velocity  as  the 
"psrent"  media.  This  is  justified  by  the  fact  that  within 
the  flow  field  our  one-dimensional  assumption  does  not 
account  for  variations  in  velocity  across  the  cross- 
section. 

c.  During  the  burning  process  the  gas  phase  gains  momentum 
equal  to  oiUp  while  the  particle  phase  loses  the  same  amount. 

Having  made  these  assumptions  we  will  now  derive  the  momentum  equa- 
tion for  the  gas  phase  which  equates  the  time  rate  of  increase  of 
momentum  in  a control  volume,  to  the  forces  acting  on  the  control 
volume  (positive  if  acting  in  the  positive  x direction).  These 
forces  include  the  reaction  of  the  net  momentum  flux  through  the 
main  entrance  and  exit  of  the  control  volume,  the  force  due  to  tKe 
momentum  flux  associated  with  mass  addition  from  burning,  the  r.  * 
action  due  to  momentum  loss  associated  with  mass  passing  through  the 
wall,  the  pressure  gradient  force,  the  Interaction  force  of  the 
particles  acting  on  the  gas  and  the  force  from  the  outside  wall,  or 

iiSAu).  + Mu  A - (o)  u A)  - A - fA  + F K <A5) 

3t  dx  p W w 3x  w 

where  the  term  u is  the  x component  of  velocity  of  the  gas  leaving 
the  control  volume  through  the  wall. 


Equation  (A5),  after  combining  with  the  gas  continuity  equation, 

Eq.  (A2),  can  be  written  as 

°Dt’"lx+Fw"F  + “(Vu)  “ “v(Vu)  (A6) 
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Similarly,  the  momentum  equation  for  the  particle  phase  is 
3(o  Au  ) 3(o  Au^) 

- (wu  A)  - (w  u A)  + FA  + F A (A7) 

3t  3x  p wp  wp  wp 

which  upon  simplification  becomes 
D^u 

a 77— ■ F + F “ u»  (u  ~u  ) (A8) 

p Dt  wp  wp  wp  p 


where  the  term  u„p  is  the  x component  of  the  velocity  of  the 
particles  passing  through  the  wall. 


In  Eqs.  (A6)  and  (A8),  the  interaction  force  between  the  gas 
and  the  particles,  F,  in  general,  includes  four  types  of  forces, 
namely,  the  viscous  drag  force,  the  pressure  gradient  force,  the 
apparent  mass  force  due  to  the  .acceleration  of  the  gas  surrounding 
the  particles  and  the  force  due  to  nonsteady  flow.  Note  that  the 
drag  force  does  not  include  the  pressure  gradient  effect.  Often, 
the  drag  force  determined  experimentally  contains  both  the  effects 
of  viscosity  and  pressure.  Special  care  must  be  taken  to  separate 
these  effects  in  using  the  present  set  of  equations.  A detailed 
discussion  on  these  forces  may  be  found  in  Hinze  [Al],  Rudinger  [A2], 
Pal  [A3]  and  Willis  [A4].  Let  us  follow  the  approach  of  Hinze,  and 
consider  the  force  acting  on  a spherical  particle,  given  by 
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If  we  now  define  the  number  of  particles  per  unit  volume  of  mixture, 
N,  as 


p x(Volume  of  a particle) 
P 

6c 
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where  F represents  the  interaction  force  without  the  pressure 
gradient  tern.  As  the  value  of  _iL  becomes  small,  the  apparent  mass 

PP 

force  and  the  force  due  to  nonsteady  flow  become  negligible.  It  is 
assumed  for  the  problem  presented  here  that  the  values  of  p/pp 
encountered  are  small  enough  so  that  these  terms  can  be  neglected 
leaving  us  with  an  expression  for  F which  includes  only  the 
viscous  drag  force,  namely: 


F - JL 
PP 


K]_  iu’u pl<u-u p) 


where 


3C.o 

K - --^-E 
1 AD 


Setting  u„  equal  to  u and  Uyp  equal  to  Up,  Eqs.  (A6),  (A8)  and 
(A9)  combine  to  become 

5?  + ^ H • - c l(F-V  ' “<Vu)+“w(Vu)1  <A10> 


and 
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(All) 


respectively.  Solving  eq.  (A10)  for  -r*-  and  substituting  into 
Eq.  (All)  yields:  X 
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(A12) 
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The  term  will  now  be  neglected  for  ainq>licity. 

Dt 


With  this  term 


neglected,  the  particle  equations  are  weakly  coupled  to  the  system,  i.e. 
the  particle  equations  do  not  contain  derivatives  of  the  gas  variables 
u,  o,  and  E.  The  ratio  of  this  term  over  the  first  term  of  A12  Is 
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which  is  much  less  than  one  when  p/p  and  e are  both  small.  For  Instance, 

p 

in  the  numerical  example  treated  later,  this  ratio  is  less  than  0.15  for 
most  points  In  the  flow  field.  However,  near  the  breech  and  during  the 
time  immediately  after  initiation,  this  ratio  may  be  larger  than  o.5  and  the 
omission  of  the  Du/Dt  term  may  be  a poor  approximation.  An  attempt  should 
be  made  to  retain  this  term  in  future  refinement  of  this  work.  Eq.  (A8)  is 
then  reduced  to 
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Equations  (A10)  and  (A13)  are  the  momentum  equations  used  in  the  computer 
code  for  the  gas  and  particle  phases  respectively. 


In  the  derivation  leading  to  Eq.  (A10) , we  have  implied  that  the 
pressure  force  acting  on  the  particles  is  - e3p /3x,  which  is  different  from 
Pal's  expression  of  -3(pe)/3x.  We  feel  that  our  approach  is  more  realistic 
for  the  present  problem.  In  order  to  get  a better  feel  for  this  term  let  us 
examine  in  more  detail  the  pressure  force  acting  on  the  particles  and 
compare  our  approach  to  Pal's  [A3].  Let  us  begin  by  discussing  Pal's  treat- 
ment of  the  pressure  Itself.  In  his  approach,  he  assumes  that  the  particle 
density  pp  is  constant  (the  particle  concentration  cp  is  variable)  and  that 
the  pressure  of  the  mixture  p,  which  he  calls  total  pressure,  is  the  true 
pressure  of  the  gas.  He  defines  the  gas  partial  pressure,  p>,  as  the 
pressure  of  the  gas  of  fixed  mass  and  fixed  temperature,  if  it  were  to 
occupy  the  entire  volume  of  the  mixture.  If  the  equation  of  state  of  the  gas  is 


p ■ p(p,E)  (A14) 

then  Pg  is  defined  by, 

Pg  - (o,E)  - p[p(l-c),  E]  (A15) 

In  addition,  if  the  pressure  is  a linear  function  of  p.(p  ■ pf(E)) 
which  is  the  case  for  an  ideal  gas,  then  pg  is  given  by 

Pg  • (1-c)  pf(E)  - (l-e)p  (A16) 

The  total  pressure  of  the  mixture  is  considered  as  the  sum  of  the 
partial  pressure  of  the  gas,  pg,  and  that  of  the  pseudo-fluid  of 
particles,  pp  or 

p " pg  + PP 
From  the  above,  it  follows  that 

Pp  " eP 
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Note  that  the  definition  of  the  partial  pressure  of  the  gas,  Pg 
follows  Dalton's  law  of  partial  pressure,  l.e.  the  partial  pressure* 
of  a component  in  a mixture  is  the  pressure  that  the  component  would 
exert  if  it  were  alone  in  the  container.  However,  in  reality, 
partial  pressure  of  the  particle  phase  pp,  does  not  obey  Dalton's 
law,  l.e.  if  the  particles  were  alone  in  the  container  they  would 
not  produce  a pressure. 

Now,  we  will  proceed  to  discuss  the  pressure  force  acting  on 
each  phase.  In  Pal's  approach,  the  particles  are  considered  to  be 
smeared  and  occupy  the  complete  volume  of  the  mixture,  with  only 
their  partial  pressure,  pp,  acting  on  them.  In  our  approach,  we 
also  consider  the  particles  to  be  smeared;  however,  we  assume  that 
they  occupy  only  a portion  of  the  volume  having  an  effective  cross- 
sectional  area  which  is  being  acted  upon  by  the  total  pressure. 

The  particle  equation  of  motion  as  derived  by  Pal  is  of  the 
form,  (the  following  discussion  also  holds  for  the  gas) 


Du  dp 

o„  -£r^-  + — , where  the  right  hand  side  contains  no  derivatives. 

The  pressure  gradient  term  In  this  equation  can  be  written  as  j 

<*X 

which  is  different  from  that  in  our  equation  of  motion,  e r*- . 

9e  dX 

In  other  words,  we  do  not  include  the  term  p — in  our  equation  of 

motion.  We  can  justify  this  by  analyzing  the  following  models. 

However,  first,  let  us  review  the  derivation  of  the  momentum  equation 
for  a single  phase  flow  with  area  change. 


C=  > + 
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Figure  A1  - Control  Volume  of  a Single  Phase  Flow 

Considering  the  control  volume  shown  in  Fig.  (Al),  the  net  pressure 
force,  pn,  acting  on  this  control  volume  is  in  the  x direction  and 
given  by 
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Pn  - + p A + (p  + f£  y)  dx-(p  + If  dx)<A  + dx) 


where  the  second  term  on  the  right  hand  side  is  the  pressure  force 
on  the  side  wall.  After  neglecting  terms  of  second  order  in  dx, 
the  above  equation  becomes 


Pn  “ - A 


iz 

3x 


dx. 


Note  that  there  are  no  first  order  contributions  from  the  change  in 
area  dA  to  the  net  pressure  force;  the  first  order  effect  of  the 


pressure  force  due  to  ~ on  the  right-hand  surface  is  cancelled  by 
the  pressure  force  on  the  side  wall. 

For  the  particles  in  a two-phase  flow  we  have  a similar 
situation  (Fig.  A2).  Let  A be  the  cross-sectional  area  of  the 
mixture,  A be  the  equivalent  cross-section  area  of  the  particles, 
and  C2  dx  be  the  equivalent  length  of  the  particles  such  that 
ei  E2  - e. 
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Figure  A2  - Control  Volume  for  a Particle  - 
Gas  TVo-Phase  Flow 

The  net  pressure  force  on  the  particles  is  then 

. e2dx  3(elA>  So  r 

p„  - pEiA  + <p  + ■&—' 1 — 3T'ix  - <p  + £ dx) 

which  after  neglecting  second  order  effects  becomes 

p ■ - 1^-  A e dx  . 
rn  3x 

Note  that  in  our  formulation,  we  treat  the  change  in  e as  a change 
in  area,  resulting  in  a net  pressure  force  on  the  particles  that  does 

not  Include  a term  of  the  form  1^. 
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The  same  conclusion  may  also  be  reached  by  considering  an  extreme 
case,  where  all  of  the  particles  are  concentrated  in  a "rod"  of  area  A e 
and  length  dx,  as  shown  in  Fig.  A3. 


Figure  A3  - Particle  Rod  Model 


The  net  pressure  force  acting  on  the  ‘'particle  rod”,  pnf,  is  then 

PAt  + (p  + U d*  - (p  + l^x)tAE  + 1^,‘1 


nr 


or  neglecting  second  order  terms 
p _ ■ - Ae  dx 


nr 


dx 


Note  that  the  pressure  acting  on  the  particle  phase  in  this  model  acts 
on  an  area  eA,  while  the  pressure  acting  on  the  gas  phase  acts  on  an  area 
(l-c)A.  Thus,  for  the  rod  model,  the  pressure  force  acting  on  the 
particle  and  gas  phases  are  respectively  - Ac  i£  and  _ wi_e\lL 

3x  v 'dx’ 

Recognizing  that  in  this  model_the  Interaction  force  does  not  contain  a 
pressure  gradient  term  (F  not  F must  be  used)  we  have  consistency  between 
the  simple  rod  model  and  the  complicated  continuum  model. 

3.  Energy  Equations 


Before  deriving  the  energy  equations  for  the  gas  and  particle  phases, 
let  us  review  the  general  energy  equation  for  a control  volume.  The 
first  law  of  thermodynamics  for  a control  volume  which  is  fixed  in  space 
(Eulerian  approach)  is 
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where 


^shaft 
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rate  of  heat  added  to  the  control  volume  (energy  per 
unit  time) 

power  delivered  from  medium  inside  the  control  volume  to 
the  outside,  due  to  rotating  mechanical  members  (work 
per  unit  time) 

power  transferred  from  medium  within  the  control  volume  to 
adjacent  outside  medium  through  shear  force. 

Note  that  the  outside  medium  must  be  in  motion  to  have  power 
delivered;  if  the  outside  medium  is  stationary,  the  power 
will  be  zero  even  though  there  is  force  transmitted. 


f [ ]dV  ■ the  time  rate  of  increase  of  the  internal  and  kinetic 
energy  of  the  medium  inside  the  control  volume. 


f [ ]dA  - the  net  enthalpy  flux  and  kinetic  energy  flux  going  out 
of  the  control  surface,  (work  per  unit  time). 

E ■ internal  energy  of  gas  per  unit  mass  of  gas 

H * E + ^ ■ enthalpy  of  gas  per  unit  mass  of  gas. 


The  energy  equation  for  the  gas  phase  is  then 

2 
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(A17) 


where  the  total  dissipation  $ has  been  broken  down  into  a boundary  dissipation 
term  and  an  internal  dissipation  term,  (see  Appendix  B)  namely: 


where 


*B  " ^“we  “ »>  “ F<«p  _ u> 


The  term  E is  the  internal  energy  of  the  combustion  products  calculated 
at  the  flame  temperature  and  is  not  included  in  q.  Care  must  be  taken 
in  supplying  the  value  of  q to  insure  that  £ is  not  included.  In  writing 
Eq.  (A17),  we  have  assumed  that  the  values  of  q and  E do  not  include  the 
effects  of  the  velocity  of  the  particles  and  the  flow  work.  Resulting 
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from  the  creation  of  gas  mass,  thus  these  terms  are  treated  separately  in 
the  5th  term  on  the  right  hand  side  of  Eq.  (A17) . It  is  felt  that  the 
kinetic  energy  attributed  to  the  velocity  of  the  combustion  products  as 
they  move  away  from  the  surface  of  the  burning  propellant  is  included  in  q 
and  need  not  be  Isolated.  The  slkth  term  on  the  right  hand  side  accounts 
for  the  energy  lost  by  gas  passing  through  the  wall.  The  last  term 
represents  the  flow  work  associated  with  particles  leaving  through  the  wall. 

When  calculating  the  work  done  on  the  surroundings,  it  should  be 
noted  that  there  is  no  work  done  by  the  wall  friction  force  since  the 
wall  is  stationary;  however,  the  interaction  force  F (remembering  that 
in  accordance  with  the  simplified  model  the  interaction  force  is  F and 
not  F)  does  do  work  since  the  gas  particle  boundary  has  a velocity  equal 
to  Up.  To  clarify  this  point,  let  us  first  consider  a viscous  fluid 
flowing  next  to  a stationary  solid  wall  (see  Fig.  (A4)).  The  velocity  of  the 
fluid  in  contact  with  the  wall  must  be  zero.  The  velocity  Increases  from 
zero  at  the  wall,  through  the  boundary  layer  to  its  "free  stream"  value  of  u. 
If  our  control  volume  is  taken  with  the  solid  wall  as  a boundary,  then  the 
shear  force  F§  does  no  work  on  the  medium  outside  the  control  volume,  be- 
cause the  boundary  of  the  medium  is  at  zero  velocity.  The  viscous  stress 
in  the  boundary  layer  within  the  control  volume  will  dissipate  mechanical 
energy  into  heat,  but  this  energy  transfer  is  within  the  control  volume  and 
does  not  represent  energy  transfer  across  the  control  volume  boundary.  If 
the  boundary  of  the  control  volume  is  another  fluid,  then  the  rate  of  work 
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done  by  the  medium  to  the  outside,  p * Is  Fs  u^  Adx,  where  ujj  is  the 
velocity  of  the  common  boundary  of  the^fuids,  not  their  free  stream 
velocities.  If  the  boundary  is  a solid  moving  at  a velocity  u^,  then  the 
rate  of  work  done  is  also  Fg  u^,  Adx. 

In  our  present  case,  the  particle  phase  is  "solid",  therefore  the  rate 
of  work  done  by  the  gas  on  the  particle  is  F Up  Adx  and  conversely  the  rate 
of  work  done  by  the  particles  on  the  gas  is  -F  Up  Adx. 

If  we  now  define  Q and  Qw  as  the  total  energy  addition  terms  due  to 
burning  and  energy  transfer  through  the  wall,  respectively,  as 
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(A18) 


(A19) 


and  simplify  Eq. (A17)  by  subtracting  Eq.(A2)  multiplied  by  (E  + u / 2)  and 
Eq.(AlO)  multiplied  by  u,  we  arrive  at  the  final  version  of  the  energy 
equation  for  the  gas  phase,  namely: 


DE  + 2.  £u  m JL 

Dt  p 3x  a 3x  o 


[ + u(F  - Fw)  - Fup  + 


- wfuu  - 4 u^  + E)  + u)  (uu  -■^■u^  + E)~w  -E-  - ~ 1 (A20) 

p2  w w 2 wp  pp  A 3x  1 

For  the  current  problem,  is  assumed  to  be  negligible.  Similarly, 
the  particle  energy  equation  can  be  written  as 

°n(qn  + “ - Fu + TT  [P_eA(E  + k + ^7lP„EA  Un^En  + \ Un  + ^ 1 
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where  we  have  defined  the  enthalpy  of  the  particle  phase,  Hp,  as 

H - E + -2-  (A22) 

P P Pp 

and  the  enthalpy  of  the  particles  leaving  through  the  wall,  H , as 

p ^ 

H - E + ~ (A23) 
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If  we  now  define  Qp  and  ()  , the  total  energy  addition  terms  due  to  burning 

and  energy  transfer  with  the  wall,  respectively,  as 
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and  simplify  Eq.(A21)  utilizing  Eqs.(A4)  and  (All)  we  arrive  at  the 
final  version  of  the  energy  equation  for  the  particle  phase;  namely: 
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A simpler  form  of  the  particle  energy  equation,  (A26),  results  if  we 
assume  that  the  pressure  term  in  Eqs.(A22)  and  (A23)  may  be  dropped  and 
the  flow  work  associated  with  the  particles  phase  is  accounted  for  by  the 

term  e Au  , This  is  equivalent  to  dropping  the  pressure  contribution 
in  the  last  two  terms  of  Eq.(A21)  and  replacing  the  term  -2_  (eAu  p)  with 


cAu  |E. 
p 3x 


The  resulting  particle  energy  equation  is  then 


3x 
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(A27) 


It  should  be  noted  that  Eqs.  (A26)  and  (A27)  are  presented  here  for  the 
purpose  of  completeness  only.  They  are  uncoupled  from  the  rest  of  our 
system  and  do  not  enter  into  the  calculations. 
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Appendix  B*.  Derivation  of  One-Dimensional  Tube  Flow  Equations 

During  the  course  of  deriving  the  two-phase  flow  equations,  some  basic 
questions  were  raised  about  the  approximations  involved  in  the  one-dimensional 
representation  of  fluid  flow  in  a tube.  Most  engineering  textbooks  treat 
the  one-dimensional  flow  problem  by  applying  the  conservation  laws  to  the 
one-dimensional  tube  directly,  without  relating  the  resulting  equations  to 
the  three-dimensional  field  equations.  In  this  appendix,  we  shall  derive 
the  governing  equation  for  one-dimensional,  one-phase  tube  flow,  from  the 
general  three-dimensional  field  equations,  and  Indicate  the  approximations 
involved. 

The  three-dimensional  governing  differential  equations  are  (see,  for 
example,  Refs.  B1  and  B2.) 

the  continuity  equation 

to  + ^(p  Ui)  - 0,  (31) 

the  momentum  equation 


(B2) 


(B3) 


where is  the  component  of  force  per  unit  mass  due  to  external  sources, 

t ^ j is  the  viscous  stress  tensor,  Q is  the  heat  addition  per  unit  mass  per 

unit  time,  q is  the  heat  conduction  within  the  element,  and  4>  is  the 
dissipation  function  given  by 


*.IT  + 

9 2 Tij  ^ 3Xj  3x1// 


We  shall  first  make  the  approximation  that  u2  <<  u^,  u^  <<  u^  and  thus  the 
only  non-vanishing  component  of  velocity  is  u^.  In  a more  precise  manner, 
we  may  relax  the  requirement  on  U£  and  U3  and  only  require  that  their  integral 
across  the  area  normal  to  the  flow  direction  be  small,  namely: 


\\k  “3  “ ” \l  "1  * (B 

where  for  convenience  we  will  Interchangeably  use  the  notation  (x,y,z) 
corresponding  to  (x^.x^x^).  For  simplicity,  we  shall  drop  the  u2  and  u^ 
terns  right  from  the  beginning.  Equations  (Bl)  to  (B3)  then  reduce  to 


* Dr.  Alvars  Celmins  first  proposed  the  definitions  of  average  quantities 
used  in  this  Appendix  and  the  approximations  involved. 
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22.  + iJM.  . o 

at  ax 


(b6  J 


Du  , 3p 

p Dt  + 3x  ”pX"  3Xj 


(B7) 


DE  . p 3u  .if  au  . 3u  , 3u  1 1 fjl  ( \ 

Dt  p 3x  " ^ o |Txx  3x  Txy  3y  xz  3z  J 0 3x^  \®  ' 


where  the  subscript  on  u has  been  dropped.  Integrating  the  continuity 
equation,  (B6)over  the  cross  section  of  the  duct,  we  have 


dx  - 0 


(B9) 


We  may  interchange  the  order  of  the  x differentiation  bn  the  area 
integration  without  using  a double  integral  extension  of  Leibnitz's  rule, 
which  is  not  readily  available,  by  considering  the  area  integration  as  a 
single  integration.  Referring  to  Fig. (Bl)  it  can  be  shown  that  if 


Duct  wall 


Figure  (Bl)  Diagram  Showing  Cross  Section  of  the  Duct 


a parameter  r " r(y,z)  exists  such  that 

r(x,y,z,t)  - r(x,r,t) 
and 

A - A(r) 

where  T Is  any  flow  variable  (note  that  for  a circular  cross  section 
r la  the  radius  and  A(r)  - irr*) 


then 


\\  £r<u-£  (jrdA-rlc£  <B12> 

A A 

where  C is  the  boundary  of  the  cross  section  of  the  duct.  If  T 
vanishes  on  C,  we  then  have 


n £ r " ■ h a r “ «i3> 

A A 

Applying  Eq.(B13)  to  Eq.(B9),  assuming  that  u vanishes  on  C,  we  arrive 
et  the  continuity  equation  in  terms  of  average  properties,  namely: 

3tT<Ap)  + u)  - 0 (B14) 

where  the  average  density  p and  average  velocity  u are  given  by 


1 

" A 

(B15) 

1 

m — 

ap 

llA0U  “ 

(B16) 

defined  by  Eq. 

(B16)  is  the  "weighted"  average,  not  the 

u*, 

given  by 

_ 1 
" A 

JJ  u dA 

(B17) 

The  quantity  u*  does  not  satisfy  the  continuity  equation  exactly. 

The  momentum  equation,  (B7),  when  combined  with  (B6)  multiplied  by  u, 
can  be  written  as 


3(p»>  + .30>.ull  + 22. 

at  ax  ax 


pX  + 


(B18) 


Integrating  this  equation  over  the  cross  section,  and  applying  Eq.  (B12), 
we  obtain 


ft  ||  + £ |f 


/“2  “ + A JJ/  “ 


-pc  fHlAoxdA+ 


(B19) 


The  last  term  in  (B19)  may  be  integrated  as 


ll  ^ • ftp?  d>  * 

■ |j  IT  Txx  + IC('T«  *> 


+ t dz)  (B20) 

xy 


cm 

where  Green's  lemma  has  been  used.  The  last  tern  of  (B20)  is  the  axial 
component  of  the  shear  force  on  the  boundary  of  the  cross  section,  which 
will  be  designated  as  F . Defining  the  averages  p and  X by 


‘ 1 11/  “ 


p X dA 


(B21) 


(B22) 


Eq.  (B19)  can  be  written  as 


di  + \\  L TXX  dA  ‘ 3x  jf 


pu^  dA  + ~(Ap  u ) 


(B23) 


3 » M 2 

wh^xe  the  term  "^(Ao  u ) has  been  added  to  both  sides  of  the  equation. 

Combining  Eq.  (B23)  with  (B14)  and  neglecting  the  last  five  terms  in 
(B23),  the  one-dimensional  momentum  equation  becomes 

|i+;|i  + i|i.jt+4  (,24) 

3t  3x  p 3x  Ap  ' 





<T  . 


(B25) 


The  approximation  involved  in  deriving  the  momentum  equation  is  then 

3t 


’ > 

< t 


xx 


3x 


dA 


<< 


F I 

U I 


\\  (.. 


2 - -2 

)dA  - Ap  u 


<<  |Ap  u^| 


(B26) 


(B27) 


The  term  is  the  frictional  force  per  unit  axial  length  between  the 
fluid  and  the  wall;  it  may  be  measured  experimentally,  or  estimated. 

The  conservation  of  energy,  in  differential  equation  form,  after 
neglecting  U2  and  U3,  is  given  by  Eq.  (B8).  Adding  Eq.(B6)  multiplied  by 
E to  Eq. (B8)  multiplied  by  0 and  integrating  over  the  cross  section,  we 
arrive  at 


_3_ 

3x 


(pu)  + u 


3p  , 3u 

-*•  + T r — 

3x  xx  3x 


, 3u  , 3u.  . , 

+ T r — + T — ) dy  dz 

xy  3y  xz  3z 


(B28) 


The  last  three  terms  in  Eq. (B28)  represent  the  viscous  dissipation, 
The  last  two  terms  will  b » simplified  according  to  the  following: 


l 


(x 


*y 


3u  , 3ux.  . 

— + t — )dy  dz 

3y  xz  3z 


3(Txvu)  3fT 
S' + — 


ay 


*1 


(t  u dz  - x 
xy  xz 
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' : * ■ 

■ 


s 


or 


♦"♦B+  *i 


(B29) 


where 


B - F u - F u 
v wa  w 


III 


(u  - u) 


+ T 


XX 


dul 

3xJ 


dA 


(B30) 


and  u^a  la  the  velocity  of  the  fluid  at  the  wall. 


The  first  tern  on 


the  right  hand  side  of  Eq.(B29),  represents  the  dissipation  of  energy 

through  wall  friction  (boundary  dissipation).  The  term  $i  represents 
Internal  dissipation.  The  reason  for  separating  the  total  viscous 
dissipation,  4»,  Into  these  two  components  is  to  conform  to  the  accepted 
engineering  practice  of  utilizing  a wall  friction  coefficient  for  one- 
dimensional  problems.  It  must  be  realized  that  the  application  of  a steady 
state  friction  coefficient  to  an  unsteady  problem  may  Involve  substantial 
error. 

If  we  now  make  the  following  definition  for  average  quantities, 


_1_ 

if  pE  dA 

(B31) 

A?  , 

U 

_1_ 
Ap  , 

fj  (PQ  - ^)dA 

(B32) 

and  neglect  certain  terms,  which  involve  the  following  approximations. 


| ||  pEu  dA  - P Eu  A | 

« |p  Eu  A| 

(B33) 

j ||  pu  dA  - p uA  << 

Ip  ua| 

(B34) 

. 

I 

1 

« i '“H*i 

1 

(B35) 

Equation  (B28)  reduces  to 

llM-Il  ilAoj-al  ■ pAQ  - Ap  ~ + F (u 
at  ax  ax  w v 
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- u)  + 


(B36) 


~ ' ’ * * ? • > 

■ . ’ 

►V  V-'.  - . -dT*  . 1 “*•  •*-  * 


St 'StU&jgfeKa:  - --  - - if B 1 


Combining  (B36)  with  (B14)  and,  (B24)  multiplied  by  u,  we  obtain. 


_D 

Dt 


4 7-  (Au  p)  + X u + Q 


<J>. 


(B37) 


If  the  dissipation  function  <t>  had  not  been  separated  as  in  F,q.(B29), 
Equation  (B37)  would  then  become 


£ r,  <*=  ’> 


+ 


F 

x u + q + 4 

Ap 


(B38) 


Under  certain  situations  in  Eq.(B37)  may  be  neglected.  It  is 
realized  that  4>  must  at  all  times  be  positive  and  if  the  problem  is  such 
that  either  $3  or  are  negative  this  assumption  must  be  reexamined. 

For  all  practical  problems,  (uy  - u)  and  Fw  are  opposite  in  sign,  therefore 
the  term  $3  is  positive.  Within  the  boundary  layer  where  ax^/ay  and 
^xz/^z  are  lar8e  a°d  positive,  the  term  (u  - u)  is  also  positive,  therefore 

$1  is  also  positive. 

The  one- dimensional  equations  are,  in  general,  good  for  flows  of  close 
to  uniform  velocity  distribution.  The  viscous  force  is  then  negligible  ex- 
cept near  the  solid  wall,  where  a thin  layer  of  boundary  conditions  exist. 
This  boundary  layer  is  responsible  for  the  wall  friction  term  Fw. 

A typical  set  of  one-dimensional  tube  flow  equations  may  be  found  in 
Ref.  (B3). 
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APPENDIX  C - COMPUTER  CODE 


This  appendix  contains  a cimplete  description  of  the  computer 
code  TWOFLO  including  input  and  output  information.  A brief  explana- 
tion of  the  function  of  each  subroutine  used  in  the  code  is  presented 
along  with  instructions  for  altering  the  user  specified  routines. 

The  nomenclature  used  in  the  actual  code  is,  wherever  possible,  con- 
sistent with  this  report  limited  only  by  the  fact  that  Greek  symbols 
and  lower  case  letters  are  not  available  in  Fortran  IV  the  code 
language.  TWOFLO  is  written  to  accept  any  consistent  set  of  units  or 
can  be  run  nondimensionally . 

I.  Description  of  Fortran  Variables 

The  following  is  a list  of  the  major  FORTRAN  variables  used  in 
the  TWOFLO  code.  Those  variables  followed  by  a dash,  e.g.  UP-,  are 
suffixed  with  several  different  qualifying  symbols,  point  numbers, 
point  letters,  variables  of  differentiation,  and  so  forth.  For  an 
explanation  of  these  symbols,  see  3 below. 

1.  Flow  variables 


CG- 

at 

c 

EG- 

s 

E 

EP- 

m 

EP 

ES- 

m 

PG- 

m 

p 

RG- 

m 

p 

RP- 

m 

Pp 

TG- 

91 

i 

UG- 

■ 

u 

UP- 

X 

Up 

XX- 

m 

X 

Miscellaneous 

variables 

GC- 

m 

GC 

GM- 

m 

GM 

GE- 

m 

GE 

PC- 

m 

PC 

PM- 

m 

PM 

PE- 

m 

PE 

F- 

m 

F 

FWG- 

m 

Fw 

FWP- 

m 

Fwp 

QG- 

m 

Q 

QP- 

m 

QP 

QWG- 

m 

Qw 

QWP- 

m 

Qwp 

W- 

m 

<D 

WWG- 

m 

UJy 

WWP- 

m 

75 


DTT- 

- At 

Dt- 

= At 

DX- 

= Ax 

3.  Suffixes  - suffixes  are  used  to  describe 

a.  mesh  points,  i.e.  SP(K, I)  where 

K - time  line  (1,2) 

I * point  number 

b.  a particular  point  in  the  iteration  scheme  i.e.  SP1,  SPA, 
etc.  where  1 and  A are  point  designations. 

c.  derivatives,  i.e.  SPX-,  SPT-  and  GEUG- 

SPX-  = 3op/3x 
SPT-  - 3cp/3t 
GESP-  - 3GE/3op 

(suffixes  X and  T represent  x and  t respectively) 


II.  Description  of  subroutines  * 

1.  ABC(L,K1,K2 ,K3, 11,I2,I3,M1,M2,M3,J2,IC, ID.NNW1.NNW2, TT1 , TT2 , 
TT3,TTT,NXT) 

This  subroutine  is  designed  to  interpolate  for  the  proper- 
ties at  the  base  of  the  characteristics  when  only  the  gas  phase  is 
present.  When  quadratic  interpolation  is  required,  a Lagrange  inter- 
polating polynomial  of  order  two  is  used:  namely 

(x-x2)  (x-x3)  (x-x^  (x-x3) 

" (x1-x2)  (x-j-x-j)  yl  + (x2-xx)  (x2-x3)  y2 

(x-x1)(x-x2) 

+ (x-j-x^  (x3-x2)  y 3’ 

where  y(x)  is  the  particular  property  to  be  calculated,  and  y^,  y2 
and  y3  are  the  values  of  that  property  at  the  mesh  points. 

The  terms  in  the  call  list  are 


K1 


1 

2 

3 

4 

5 

6 


calculate  properties  at  points  A,B  and  C 
" " " point  A 

fl  tt  ff  II  g 

" " " points  A and  C 

" " " " B and  C 

" " " point  C 


* Note:  For  all  subroutines,  L ■ 1 Indicates  that  the  point  is  located 
behind  the  bullet;  L ■ 2 indicates  that  point  is  located  in 
front  of  the  bullet. 
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K2 

s 

K2 

i 

K3 

11,12,13 

Ml , M2 , M3 

3 

J2 

3 

J2 

= 

1C 

3 

ID 

3 

NNW1 

« 

NNW2 

• 

TT1, TT2, TT3 

HI 

NXT 


1 C properties  are  assigned  to  be  the  same  as  point  (M1,IC) 

1 C properties  are  interpolated 

1 C properties  are  not  calculated  (used  after  first 
iteration  when  K2  ■ 1) 

point  numbers 

line  number  of  points  11,12,13  respectively.  (1  old  time 
line 2 new  time  line  ) 

1 quadratic  interpolation  is  used  between  points  11,12  and  13 

2 linear  interpolation  is  used  between  points  11  and  12. 
(K2*1)IC  is  the  point  number  of  the  mesh  point  correspond- 
ing to  point  C. 

not  used 

tl  ft 

tl  II 


* times  associated  with  points  11,12,  and  13  respectively 
if  interpolation  is  with  respect  to  time  instead  of 
position. 

* not  used 

* 1 interpolation  along  a constant  time  line  using 

position 

■ 2 interpolation  using  time 


2.  ABCD(L,K1,K2,K3,I1,I2,I3,M1,M2,M3,J2,IC,ID,NNW1,NNW2, 

TT1,  TT2,TT3, TTT.NXT) 

This  subroutine  is  quite  similar  to  ABC  except  here  the 
properties  of  both  phases  are  calculated.  The  terms  in  the  call  list 
which  are  different  from  ABC  are 


1 

calculate  properties 

at  points 

A,B,C  and  D 

2 

II  II 

II 

II 

A,C  and  D 

3 

II  tl 

II 

II 

B,C  and  D 

4 

II  II 

II 

point 

A 

5 

II  II 

II 

II 

B 

6 

• 1 II 

II 

19 

B 

7 

II  II 

II 

If 

D 

8 

II  II 

II 

It 

C 

9 

• 1 II 

II 

II 

C 

20  do  no  calculating  just  print  A,B,C  and  D properties 

K2  - 1,4,5  properties  at  point  C correspond  to  properties  at 
point  (Ml.IC) 

2,3  properties  at  point  D correspond  to  properties  at 
point  (Ml, ID) 

ID  - the  point  number  of  the  mesh  point  corresponding  to 
point  D if  K2  ■ 2 or  3. 

3.  ADISC(L,ID,K1,K2,K3,K4,K5) 

This  routine  has  two  purposes:  it  handles  discontinuities  in 

cross-sectional  area,  A(x)  , and  is  used  also,  either  in  conjunction 
with  INRTBN  or  simply  by  itself,  to  handle  the  gas-particle  interface, 
when  that  interface  is  treated  as  a discontinuity  in  gas  concentration 
(o'*.  Both  of  these  purposes  are  accomplished  with  much  the  same 


ll 


computational  procedure.  There  are  two  main  differences!  however. 

The  first  is  that,  in  the  case  of  the  gas-particle  interface,  no 
particle  properties  are  calculated  in  front  of  the  discontinuity, 
since  no  particles  are  present  there;  the  second  is  that,  in  the 
case  of  the  gas-particle  interface,  the  u-c  characteristic  may  inter- 
sect with  the  path  of  the  bullet,  and  if  it  does,  the  properties  at 
the  base  of  that  characteristic  must  be  obtained  by  interpolation 
along  the  bullet  path. 

The  terms  in  the  call  list  are: 

ID  ■ point  number  of  the  point  on  the  left  hand  side  of  the 
area  discontinuity 

K1  - 0 normal  Area  Discontinuity 
* 1 Gas-particle  interface 
K2  s not  used 
K3  = " " 

K4  - " " 

K5  = " ” 

4.  ADT(C1,C2,C3,U1,U2,U3,X1,X2,X3,UP1,UP2,UP3,K1,K2,KK3) 

This  subroutine  calculates  the  At  associated  with  each  point. 

The  derivation  is  as  follows: 

4 

'T' 

s'  / ' 

/ ' 
t \ 


Figure  Bl.  Plot  of  characteristic  grid  used  to  calculate  At. 

Referring  to  Fig.(Bl),  we  first  calculate  At  (DTI) 
associated  with  the  (u+c)  characteristic  direction  by  writing 

V*1  <u4+c4)  + <UJ+Ci> 


Vx2  Vu2 


Rewritten,  these  equations  are 


xx  + DTI  • - (u^  + c4  + ux  + Cl) 
x2  + DTI  • j <uA  + u2) 


•M&k**.  . .Mtk.  -YHhirtlMtu  n ~rt 


Subtracting  these  two  equations,  we  obtain 

0 - (Xl-x2)  + DTI  • | (c4  + cx  + ux  - u2) 
Tf  we  now  assume  that  * c2»  we  have 


DTI  - 


2(*2  - Xj) 


=2  + C1  * “l  * u2 


A similar  argument  leads  to  At(DT3)  associated  with  the  (u-c) 
characteristic,  namely: 


DT3  ■ 


The  terms  in  the  call  list  are 


C1,C2,C3, 

U1,U2,U3 

X1.X2.X3 

UP1.UP2.UP3 


CG 

UG 

XX 

UP 


at  points  1,2  and  3 respectively- 

ft  tl 
It  II 


K1  ■ 1 interior  point 

2 interface  between  gas  and  two-phase  region 

3 area  discontinuity 

K2  - 1 smallest  value  of  At  resulting  from  u+c  and  u-c 
is  chosen 

2 At  resulting  from  u+c  is  chosen 

3 " u-c  " " 

KK3  * not  used. 


5.  AIN  ( IFILE , IA , KK1 , KK2 , KK3 , KK4 ) 

This  subroutine  is  used  to  read  in  data  from  a dump  tape. 
The  terms  in  the  call  list  are 


IFILE  - number  of  the  file  from  which  data  is  read. 
IA  ■ data  is  assigned  to  mesh  point  (IA,  ). 

KK1  ■ not  used 
2 m ••  *• 

3 - »•  »• 

A a tl  ft 


6.  AOUT  (IFILE,  IA.KK1.KK2.KK3.KK4) 

This  subroutine  is  the  same  as  AIN  except  that  it  dumps  data  on 
IFILE 
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7.  AR(XX,Z1,Z2)  USER  SPECIFIED 

This  function  subroutine  calculates  the  cross-sectional  area 
of  the  duct,  A,  given  the  x location. 


8.  ARX(XX,Zl,Z2)  USER  SPECIFIED 

This  function  subroutine  calculates  3A/8x  given  x. 

9.  BNPTG(L, II, 12, 13.MAXIT,  Jl, J2, J3, JA) 

This  subroutine  calculates  the  gas  phase  properties  at  boundary 
points  using  one  phase  equations 

The  terms  in  the  call  list  are: 

11,12,13  • base  mesh  points  (point  being  calculated  is  II 
if  J1  - 1 or  13  if  J1  ■ 2) 

MAXIT  • maximum  number  of  iterations 
J1  ■ 1 left  boundary 
2 right  boundary 

J2  ■ 1 specify  velocity  on  the  boundary 
2 " pressure 

J3  • 1 normal  boundary 

■ 2, 3, A special  boundary  points 

JA  - 1 J2-1  the  bullet  momentum  equation  is  used 

2 J2-1  the  velocity  (bullet)  is  zero. 

3 boundary  conditions  are  determined  prior  to 
entering  BNPTG(not  in  GBCOND) 

10.  BNSH(L,K1,JA, MAXIT, JJ1, JJ2, JJ3) 

This  subroutine  controls  the  calculation  of  a shock  reflecting 
from  a boundary 

The  terms  in  the  call  list  are 

K1  - shock  number 

JA  ■ 1 bullet  is  moving 

2 bullet  is  stationary 

MAXIT  « maximum  number  of  iterations 

JJ1  = not  used 

JJ2  - " " 

JJ3  = " " 
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11.  C0MEQ1(L,M,K1,K2,K3,Z1,Z2,Z3/ 

Using  the  compatibility  equation  along  the  characteristic 
direction  dx/dt  ■ Up,  Eq. (50), this  routine  solves  for  Up  at  point  4, 
given  values  of  the  other  material  properties. 

The  terms  in  the  call  list  are: 

M - 1 
K1  - 1 

Z1  ■ not  used 

Z2  - " " 


12.  C0MEQ2(L,M,K1,K2,K3,Z1,Z2,Z3) 

This  subroutine  uses  the  particle  energy  equation,  Eq. (52)  written  in 
finite-difference  form  to  solve  for  Ep  at  point  4,  given  values  of  the 
other  material  properties. 

Note:  This  routine  is  not  used,  because  the  particle  energy 
equation  is  uncoupled  from  the  rest  of  the  system. 

13.  COMEQ3(L,M,Kl,K2,K3,N2,Z2,Z3) 

This  subprogram  calculates  op  at  point  4 given  values  of  the 
other  material  properties  at  point  4.  The  particle  continuity  equation 
Eq.(51),  written  in  finite-difference  form,  is  used. 

The  terms  in  the  call  list  are: 

M - 1 

K1  - 1 

K2  - not  used 

K3  - J2  in  GNPT 

N2  ■ number  of  point  being  calculated  (point  4) 

Z1  - pot  used 

z2  „ " » 

14.  C0MEQ4  (L,  M,  K1 , K2 , K3 , Zl,  Z2 , Z3) 

This  routine  uses  the  compatibility  equation  along  the  direction 
dx/dt  - u,  Eq.  ( 47) . 

The  terms  in  the  call  list  are: 

M • 1 two  phase  compatability  equations  is  used 

2,3  gas  phase  " " " " 

K1  ■ 1 equation  (47)  is  solved  for  EG4 

2 " (47)  " " " " 

and  [3E/3o]4  and  [3E/9u)A  (see  NEWTON) 
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K2  ■ not  used 

K3  - 1 PM4 , PC4 , GCA , GE4  are  recalculated 

2 " " " " are  not  recalculated 

Z1  - not  used 

Z2  - " " 

Z3  - " " 


15.  C0MEQ5(L,M,K1,K2,K3, Zl, Z2,Z3) 

This  subroutine  solves  the  compatibility  equation  along  the 
direction  dx/dt  - u + c,  Eq.(48). 

The  terms  in  the  call  list  are: 


M 

M 

K1 


1.3 

2 

1 

2 


K2 

K3  - 
Zl  - 
Z2  - 
Z3  ■ 


1 

2 


two  phase  compatability  equation  is  used 
gas  " " " " » 

equation  is  solved  for  either  SG4  or  UG4 
(see  K2) 

equation  is  put  into  the  form  g(o,u)  - 0 
and  the  derivatives  3g /3a  and  3g/3u  are 
calculated  (see  NEWTON) 
solve  for  SG4(M»1,3)  or  RG4(l^2) 
solve  for  UG4 

S.A.  C0MEQ4 
not  used 

II  It 


16.  C0MEQ6(L,M,K1,K2,K3,Z1,Z2,Z3) 

This  subroutine  solves  the  compatibility  equation  along  the 
direction  dx/dt  ■ u - c,  Eq.(49). 

The  terms  in  the  call  list  are: 


M - 1,2 

- 3 


K1 


K2 


1 

2 

1 

2 


K3  - 

Zl  - 
Z2  - 
Z3  - 


two  phase  compatability  equation  is  used 
gas  phase  " ,f  " " 

equation  is  solved  for  either  SG4,  RG4  or 
UG4  (see  K2) 

equation  Is  put  in  the  form  f(a,u)  * 0 and 
the  derivatives  3f/3o  and  3f/3u  are  calculated, 
solve  for  SG4(M-1,2)  or  RG4(M*3) 
solve  for  UG4 

S.A.  C0MEQ4 

not  used 

II  It 


17.  C0MEQ7 (L , M , K1 , K2 , Zl , Z2 , Z3) 

This  routine  calculates  the  regression  distance  (Z). 
L is  the  only  term  in  the  call  list  that  is  used 
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18.  C0MEQ8(L,M,K1,K2,K3,Z1,Z2,Z3) 

This  subroutine  calculates  field  properties  when  UG4  and  UP  4 
are  zero  and  all  base  point  properties  are  identical. 

L is  the  only  term  in  the  call  list  that  is  used. 


19.  CSQ (L , M, K1 , K2 , RH , E , C , P , PE , PR , T , XA, XB) 

This  subroutine  calculates  the  sound  speed  of  the  gas 
The  terms  in  the  call  list  are: 


M 

M 

K1 

K2 

RH 

E 

C 

P 

PE 

PR 

T 

XA 

XB 


P 

E 

c 

P 

ap/3E 

3p/3p 

T 


RH  and  E must  be  input 
P,  RH,  PE,  PR  must  be  input 


not  used 

tl  It 


not  used 

ft  It 


20.  DIFF(L,N,M,K2,K3,Z1,Z2,Z3) 

This  subroutine  uses  forward,  backward  and  central  difference 
schemes  to  calculate  the  terms  3up/3x  and  3ap/3x.  Then  the  particle 
continuity  and  momentum  equations  are  used  to  calculate  the  terms 
3ap/3t  and  3up/3t  respectively. 

The  terms  in  the  call  list  are: 


N - point  number 
M - line  number 
K2  » 1 central  difference  is  used 
■ 2 backward  difference  is  used 
- 3 forward  " " " 

• 4 x derivatives  are  calculated  externally  thus  only  t 
derivatives  are  calculated. 


K3 

- M 

Z1 

■ 

not 

used 

Z2 

m 

not 

used 

Z3 

m 

not 

used 

83 


. 

■ 


21.  DIMEN (L,K1,K2,K3,K4,K5) 

This  subroutine  nondimensionalizes  all  variables  used  in  the 
program  using  B(l,20),  B(l,21),  B(l,22)  as  nondimens ionalizing  factors. 

The  terms  in  the  call  are: 


K1  - 1 nondimens lonalize  B properties 

2 nondimensionallze  points  K4  to  K5  on  the  K3 
time  line. 

3 dimens lonalize 


K2  - 
K3  - 
K4  - 
K5  - 


not  used 

time  line  number 

first  point  to  be  treated 

last  " " " " 


22.  DROP(IDR) 


This  subroutine  drops  the  IDR  point  from  the  1 time  line 

23.  DRVX(L,N1,N2,N3,K1,K2,K3,Z1,Z2,Z3) 

This  subroutine  contains  both  2 and  3 point  forward  and  back- 
ward and  central  difference  schemes. 


The  terms  in  the  call  list  are: 

N1,N2,N3  ■ points  to  be  used  in  differencing  scheme 

K1  - 1 calculates  derivatives  of  oB 
2 " " " up 

K2  ■ 1 forward  difference 

2 central  " 

3 backward  " 

K3  ■ line  number 
Z1  m not  used 

79  m 11  M 

Z3  - " " 

24.  DTQ(L, ID1, ID2, ID3, DTT) 

This  subroutine  calculates  At  for  each  point  on  the  time  line. 
The  smallest  value  after  being  reduced  by  a factor  (0.9)  becomes  the 
At  for  the  next  time  line  calculation 


The  terms  in  the  call  list  are: 


ID1 

m 

not 

used 

ID2 

m 

it 

II 

ID3 

m 

ti 

II 

DTT 

m 

At 
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25.  DZDT  ( SG , SP , EG , EP , UG , UP , PG , PP , RG , RP , TG , TP , ES , XX , DT , TT , 
K1,K2,K3,Z1,Z2,L)  USER  SPECIFIED 

This  function  subroutine  is  used  to  calculate  the  regression 
rate,  D^Z/Dt. 

The  only  requirement  in  the  call  list  is  that  the  variables 
used  in  the  equation  for  the  regression  rate  be  available. 

26.  ERH  (L , IQ , El , PI , RHI , PMEES , PMRS , XA,  TI ) 

This  subroutine  uses  the  equation  of  state  to  calculate  p 
or  E.  The  Newton  Raphson  technique  is  used  except  for  an  Ideal  gas. 

The  terms  in  the  call  list  are: 

IQ  ■ 1 calculate  p given  E and  p 
2 " E " p " p 

El  - E 
PI  - p 
RHI  - p 
PMEES  - 3p/3E 
PMRS  - 3p/3p 
XA  " not  used 
TI  - T 

27.  FQ(SG,SP,EG,EP,.UG,UP,FG,PP,RG,RP,TG,TP,ES,XX,DT,TT 
K1,K2,K3,Z1,Z2,L)  USER  SPECIFIED 

This  subroutine  calculates  the  drag  force,  F,  between  the 
particles  and  the  gas. 

The  only  requirement  in  the  call  list  is  that  the  variables 
used  in  the  equation  for  the  drag  force  be  available. 

28.  FRET(L,K1,K2,FRW,FTW,FEW,TW,RW,EW,XA,XB)  USER  SPECIFIED 

This  subroutine  computes  the  derivatives  of  f(p,E,T)  with 
respect  to  p,  E and  T. 

The  terms  in  the  call  list  are: 

K1  ■ not  used 

K2  - 1 Virial  equation  of  state. 

■ 2 Van  der  Waal  equation  of  state. 

K3  “ not  used 
FRW  - 3f/3p 

FTW  - 3f/3T 

FEW  - 3f/3E 

TW  - T 

RW  - p 

EW  - E 

XA  ■ not  used 

XB  * not  used 
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29.  FWGQ  ( SG , SP , EG , EP , UG , UP , PG , PP , RG , RP , TG , TP , ES , XX , DT , TT 
K1,K2,K3,Z1,Z2,L)  USER  SPECIFIED 

This  subroutine  calculates  the  wall  friction  force  acting  on 
the  gas,  Fwg. 

The  only  requirement  in  the  call  list  other  than  K1  is  that 
the  variables  used  in  the  equation  for  the  wall  friction  force, 

Fwg,  be  available. 

K1  ■ 1 calculate  FWg 

■ 2 calculate  FWg,  3FWg/3u  and  3Fwg/3p 

30.  FWPQ  ( SG , SP , EG , EP , UG , UP , PG , PP , RG , RP , TG , TP , ES , XX » DT , TT 
K1,K2,K3,Z1,Z2,L)  USER  SPECIFIED 

This  subroutine  is  the  same  as  FWGQ  except  that  it  calculates 
the  wall  friction  force  acting  on  the  particles,  F . 

31.  GBC0ND(L,M1,M2,K3,K4,XA,XB,PGG,UGG) 

The  subroutine  is  used  to  supply  gas  boundary  conditions. 

The  terms  in  the  call  list  are: 

Ml  ■ 1 left  boundary 
2 right  boundary 

M2  ■ 1 velocity  is  specified  on  the  boundary 

2 ii  ii  n n it 

pressure 

K3  ■ point  nunfcer  of  the  boundary  point 
K4  - 2 bullet  has  not  moved  and  u ■ 0. 

XA  * not  used 

XB  - not  used 

PGG  - pressure 

UGG  - velocity 

32.  GC(SG,SP,UG,UP,EG,EP,RG,RP,CG,PG,PP,TG,TP,ES,F,  FVG , FWP 
W,WWG,WWP,QG,QP,QWG,QWP,AR,ARX,Z1,Z2,L) 

This  subroutine  calculates  the  right  hand  side  of  the  gas 
continuity  equation,  GC. 

The  only  requirement  in  the  call  list  other  than  Z1  is  that  the 
variables  used  in  the  equation  for  GC  be  available,  (the  variables 

F, FWG QWP  are  values  generated  by  subroutines  FQ.FWGQ — 

QWPQ  respectively  and  AR  and  ARX  are  A and  3A/3x  respectively) 

Z1  * 1 calculate  GC 

- 2 " GC,  3GC/3u  and  3GC/3o 


33.  GE  (SG,  SP , UG , UP , EG , EP , RG , RP , CG , PG , PP , TG , TP , ES « f , FWG , FWP 
W , WWG , WWP , Q G , Q P , QWG , Q WP , AR , ARX , Z 1 , Z 2 , L ) 

This  subroutine  is  the  same  as  GC  except  that  the  right  hand 
side  of  the  gas  energy  equation,  GE,  is  calculated. 





34.  GM(SG,SP,  UG,UP,EG,EP,RG,RP,  CG,PG,PP,TG,TP,ES  ,F, FWG, FWP, 
W,WWG,WWP,QG,QP,QWG,QWP,AR,ARX,Z1,Z2,L) 

Phis  subroutine  is  the  same  as  GC  except  that  the  right  hand 
side  of  the  gas  momentum  equation,  GM,  is  calculated. 

35.  GNPT(L,I1,I2,I3,MAXIT,J1,J2,J3,J4) 

This  subroutine  calculates  the  properties  of  mesh  points  in 
the  two-phase  region. 

The  terms  in  the  call  list  are: 

11,12,13  * point  numbers 

MAXIT  - maximum  number  of  iterations 

J1  ■ 1 base  properties  at  points  A,B,C  and  D 
are  calculated  internally. 

■ 2 base  properties  at  points  A,C  and  D are 

calculated  Internally.  The  B properties 
are  calculated  externally. 

- 3 base  properties  at  points  B,C  and  D are 

calculated  internally.  The  A properties 
are  calculated  externally. 

- 4 base  properties  at  points  A,B,C  and  D 

are  calculated  externally. 

J2  ■ 1 normal  (Interior)  point  is  calculated. 

■ 2 point  4 is  located  using  Up.  The  gas 

compatabillty  equation  is  used  along  dx/dt 
■ u + c. 

• 3 point  4 is  located  using  up.  The  gas 
compatabillty  equation  is  used  along 
dx/dt  ■ u - c. 

- 4 two  phase  left  boundary  is  calculated. 

- 5 two  phase  right  boundary  is  calculated. 

■ 6 point  in  front  of  right  traveling  shock 

is  calculated. 

■ 7 point  in  front  of  left  traveling  shock 

is  calculated. 

J3  ■ 1 normal  general  point 

2 initial  calculation  where  interface  and 

bullet  boundary  are  calculated  simultaneously. 

■ 3 area  discontinuity  point. 

J4  • not  used. 
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36.  GNPTG(L,I1,K2,K3,MAXIT,J1,J2,J3,J4) 

This  subroutine  calculates  the  properties  of  mesh  points  In 
the  gas  only  region. 

The  terms  In  the  call  list  are: 

11,12,13  ■ point  numbers 
MAXIT  ■ maximum  number  of  iterations 
J1  ■ not  used 

J2  * 1 normal  (interior)  point  is  calculated. 

■ 2, 3, 4, 5 not  used 

• 6 point  in  front  of  right  traveling  shock  is 
calculated. 

- 7 point  in  front  of  left  traveling  shock  is 
calculated. 

J3  * not  used 
J4  - " " 


37.  GUESS (L,K1, 12, J2, J3,M, JJ5, JJ6) 

This  subroutine  supplies  the  first  guess  for  the  iterations  in 
BNPTG.GNPT  and  GNPTG. 


The  terms  in  the  call  list  are: 


i 


K1 

12 

J2 


J3 


1 first  guess  for  BNPTG  and  GNPTG 

2 " " " GNPT 

point  number  of  properties  used  as  a 
first  guess. 

1.3  (Kl-2)  first  guess  is  the  point  (M,12) 

2.4  10  (Kl-2)  first  guess  is  solution  of 

linearised  two  phase  equations 
6,7  (Kl-1)  first  guess  is  the  point  (M,I2) 

1,2, 3,4,  (Kl-1)  first  guess  is  a solution  of 

5,8,9,10  linearised  gas  phase  equations. 

not  used 


M - 

JJ5  - 
JJ6  - 


line  number  of  properties  used  as  a 
first  guess, 
not  used, 
not  used. 


38.  INDISC(L.ID) 

This  subroutine  calculates  the  complicated  singularity  that  occurs 
when  an  area  discontinuity  exists  at  the  initial  bullet  location. 


The  terms  in  the  call  list  are: 

ID  ■ point  number  of  the  left  point  at  the  area  discontinuity. 
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39.  INIT(L,IP,XB,XC,XD) 

This  subroutine  reads  in  and  prints  out  all  initial  data. 

The  terms  in  the  call  list  are: 

IP  ■ 0 read  in  and  print  out  initial  data 
J 1 only  print  out  B(  ) array. 

2 only  print  out  initial  time  line 

XB  “ not  used 

XC  - " " 

XD  ■ " " 


40.  INRTBN(L,M,K1,K2,K3, 11,12, 13, Jl, J2, J3.MAXIT) 

This  subroutine  calculates  the  gas-particle  interface  and  the 
bullet  boundary  simultaneously.  It  is  used  until  a mesh  point  is 
inserted  between  the  two. 

None  of  the  terms  in  the  call  list  are  currently  used. 


41.  INTPRN(L,TU,TL,TP,ZZ1,ZZ2,ZZ3,ZZ4) 

This  subroutine  interpolates  for  and  prints  out  the  properties 
at  a time  line  lying  between  two  calculated  time  lines. 

The  terms  in  the  call  list  are: 

TU  - time  of  upper  time  line 
TL  - " " lower  " " 

TP  ■ " " Interpolated  time  line 

ZZ1  - not  used 
ZZ2  ■ " " 

ZZ3  - " " 

ZZ4  - " " 


42.  INTSEC (XI , U1 , C 1 , T1 , X2 , U2 , C 2 , T2 , X4 , U4 , C4 , T4 , N , XB , UB , CB , TB ) 


This  subroutine  calculates  the  x and  t location  of  the  point  of 
intersection  B,  of  the  line  connecting  points  1 and  2 and  a 
characteristic  emanating  from  point  4. 

The  terms  in  the  call  list  are: 


X 

U 

C 

T 

N 


x location  of  point  _ 
gas  velocity  of  point__ 
sound  speed  of  point_ 
time  of  point_ 


1 

2 

3 


characteristic  emanating  from  point  4 is  dx/dt 

» " " " " " dx/dt 

" " " " " " dx/dt 


u + c 
u * c 
u(u>0) 


" . " M " dx/dt  - u(u<0) 
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43.  LIMPNO(Y) 

This  subroutine  drops  mesh  points,  one-by-one,  until  the  number 
of  points  remaining  is  equal  to  MAXPT  or  MAXPTB.  Each  mesh  point 
dropped  by  LIMPNO  lies  in  the  region  of  smallest  change  in  "Y" 

(with  respect  to  x)  where  Y is  specified  in  the  call  statement  and 
can  be  any  one  of  the  flow  variables.  This  helps  to  Insure  that  a 
comparatively  fine  mesh  exists  in  regions  of  large  change,  and 
that  a somewhat  larger  mesh  exists  in  regions  of  small  change. 

4A.  LINAB(L,K1,K2,K3,I1,I2,J1,J2,XP,TP) 

This  subroutine  uses  linear  Interpolation  to  calculate 
gas  properties  at  points  A or  B. 

The  terms  in  the  call  list  are: 

K1  • 1 calculate  A properties 
2 " B " 

K2  - 1 input  coordinates  of  end  points  (II, Jl)  and  (12, J2) 

■ 2 input  end  point  properties  through  common 
blocks  LIN1.LIN2.LIN3 


K3  - 1 

Interpolate  using  position 

2 

" " time,  t. 

11,12  - 

line  numbers 

Jl, J2  - 

point  numbers 

XP  - 

x location  of  point  A or  B 

TP  * 

time  of  point  A or  B 

A5.  LOCABC (Kl, II, 12,13) 

This  subroutine  calculates  the  x location  of  points 
A,B,C  and  D. 

The  terms  in  the  call  list  are: 

Kl  - 1 point  C corresponds  to  point  (1,12) 

— 2 " d " •»  •»  " 

11  ■ not  used 

12  ■ point  number 

13  ■ not  used 
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46.  NEWTON (L,M,K1,K2,K3,K4,Z1,Z2,Z3,Z4) 

This  subroutine  uses  Newton  Raphson  techniques  to  solve  com- 
patibility equations  along  characteiistics  dx/dt  = u,  u + c 

The  terms  in  the  call  list  are: 

M ■ 1 two  phase  compatibility  equations  are  used 
■ 2 gas  " " equation  is  used  alone 

(dx/dt  ■ u + c 

3 " " " " is  used  alone 

(dx/dt  * u - c 

K1  ■ not  used 

K2  * not  used 

K3  ■ 1 GC,GM,GE,PC  and  PM  are  held  constant  during  iteration 
2 GC,GM,GE,PC  and  PM  vary  during  iteration 
K4  - J2  in  GNPT-used  to  determine  the  type  of  point  being 
calculated 

Z1  - not  used 
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47.  NTROPY (L, IQ,RG, EG, SS, JJ1, JJ2 , JJ3)  USER  SPECIFIED 

This  subroutine  is  user  specified.  It's  function  is  to 
calculate  either  the  entropy,  SS,  density,  RG,  or  specific  internal 
energy,  EG  of  the  gas  phase  given  the  other  two. 

The  terms  in  the  call  list  are: 

IQ  - 1 calculate  entropy 

■ 2 " density 

• 3 " specific  internal  energy 

SS  - entropy 
JJ1  ■ not  used 

JJ2  - " " 

JJ3  - " " 

48.  PC(SG,SP,UG,UP,EG,EP,RG,RP,CG,PG,PP,TG,TP,ES,F,FWG,FWP, 

W , WWG , WWP , QG , QP , Q WG , QWP , AR , ARX , Z 1 , Z2 , L ) 

This  subroutine  calculates  the  right  hand  side  of  the 
particle  continuity  equation, PC. 

The  only  requirement  on  the  terms  in  the  call  list  is  that 
all  variables  used  in  PC  be  available. 


I 

i 
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49.  F JV (NLINE , K1 , P4 , PA, PB , PDXA , PD, XB , PDTA, PDTB) 

This  subroutine  calculates  a , o , u or  u at  point  4. 

P*x*  p » z p,x  P,t 

P4,  PA  and  PB  are  the  values  at  points  4,  A and  B of  the  variable  to 
be  differentiated.  PDXA  PDXB  PDTA  aqd  PDTB  are  the  values  of  the 
derivatives  of  the  said  variable  at  Points  A and  B. 

The  terms  in  the  call  list  are: 

NLINE  ■ not  used 

K1  - 1 calculate  x derivatives 


50.  PE(SG,SP,UG,UP,EG,EP,RG,RP,CG,PG,PP,TG,TP,ES,F,FWG,FWP, 
W,WWG,WWP,QG,QP,QWG,QWP,AR,ARX,Z1,Z2,L) 

This  subroutine  calculates  the  right  hand  side  of  the  particle 
enegry  equation,  PE. 

The  only  requirement  on  the  terms  in  the  call  list  is  that 
all  variables  used  in  PE  be  available. 

51.  PLOTTO(IT) 

This  subroutine  is  used  for  plotting  purposes  only. 


52.  PM(SG, SP,UGtUP,EG,EP,RG,RP,CG,PG,PP, TG, TP,ES,F,FWG,FWP, 

W,  WWG , WWP , QG , QP , QWG , QWP , AR,  ARX,  Z1 , Z2  , L) 

This  subroutine  calculates  the  right  hand  side  of  the  particle 
momentum  equation,  PM. 

The  only  requirement  on  the  terms  in  the  call  list  is  that 
all  variables  used  in  PM  be  available. 


53.  PPLOT(X,Y,NP,NPLOT) 

This  subroutine  is  used  for  plotting  purposes  only. 


54.  PRINTO (L , K1 , K2 , K3 , K4 , K5 , KK6 , KK7 , KK8 , KK9 , KK10) 

This  subroutine  controls  the  printout  of  calculated  data. 


The  terms  in  the  call  list  are: 


K1  - number  of  first  point  to  be  printed 
K2  - " " last  " " " " 

K3  ■ increment  by  vhich  points  are  printed  (usually  1) 

K4  ■ number  of  line  to  be  printed  (either  1 or  2) 

K5  ■ 1 print  gas  properties  only 
■ 2 print  two  phase  properties 
KK6  • not  used 
KK7  - " " 

KK8  - " " 

KK9  - " " 

KK10  - " " 
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55.  PTARNG(MSR,L) 

The  subroutine  controls  the  adding  ana  dropping  of  points 
to  maintain  a specified  time  increment. 

The  terms  in  the  call  list  are: 

MSR  - 1 always 


56.  PTQ(L,K1,K2,K3,RT,TT,PT,PTT,PTR,XA,XB)  USFR  SPECIFIED 
This  subroutine  calculates  p,  3p/3T  and  3p/3p 


The  terms  in  the  call  list  are: 


K1  not  used 

K2  ■ 1 Virial  equation  of  state 

2 Van  der  W 1 equation  of  state 

3 Ideal  gas  " " 

K3  * not  used 

Rt  - P 
TT  - T 
PT  - p 

PTT  - 3p/3T 

PTR  - 3p/3p 

XA  - not  used 

Xb  ■ not  used 

57.  QET(L,K1,K2.RT,TT,ET,XA,XB,XC)  USER  SPECIFIED 

This  subroutine  uses  the  equation  f(p,E,T)  - 0 to  calculate 
T given  p and  E or  to  calculate  E given  p and  T 

The  terms  in  the  call  list  are: 


K1 

K2 

RT 

TT 

ET 

XA 

XB 

XC 


1 calculate  T 

2 " p 

1 Virial  equation  of  state  is  used 

2 Van  der  Waal  equation  of  state  is  used 
P 

T 

E 

not  used 

II  If 
II  It 


58.  QGQ (SG , SP , EG, EP , UG , UP , PG , PP , RG , RP , TG , TP , ES , XX, DB , TT, 
Kl,K2tK3,Z,L) 

This  subroutine  calculates  the  rate  at  which  total  energy  is 
released  during  burning  of  the  propellant  Q. 

The  only  requirement  on  the  terms  in  the  call  list  is  that 
all  variables  used  in  Q be  available. 


93 


tMiUMWM— I. 


59.  QPQ(SG,SP,EG,EP,UG,UP,PG,PP,RG,RP,TG,TP,ES,XX,DT,TT, 
K1,K2,K3,ZZ,Z2,L) 

This  subroutine  calculates  the  rate  at  which  the  particles 
lose  energy  during  burning  Qp. 

The  only  requirement  on  the  terms  in  the  call  list  is  that 

all  variables  used  in  Q be  available. 

P 

60.  QUAD(Y1,Y2,Y3,XP) 

This  subroutine  uses  quadratic  interpolations  to  calculate 
the  Y property  at  location  XP  given  Y at  points  XXI, XX2  and  XX3 
(i.e.,  Y1,Y2,Y3  respectively) 

61.  QWGQ(SG,SP,EG,EP,UG,UP,PG,PP,RG,RP,TG,TP,ES,XX,DT,TT, 
K1,K2,K3,CG,Z3,CG,Z2,L) 

This  subroutine  calculates  the  total  energy  lost  by  the  gas 
phase  due  to  gas  passing  through  the  wall,  Qy  . 

The  only  requirements  on  the  terms  in  the  call  list  is  that 
all  variables  used  in  be  available. 

62.  QWPQ  ( SG , S P , EG , EP , UG , UP , PG , PP , RG , RP , TG , TP , ES , XX , DT , TT , 
K1,K2,K3,Z1,Z2,L) 

The  subroutine  calculates  the  total  energy  lost  by  the  particle 
phase  due  to  particles  passing  through  the  wall,  Q^. 

The  only  requirements  on  the  terms  in  the  call  list  is  that 
all  variables  used  in  be  available. 

63.  RHSGAS(L,K1,K2,RG,EGIPG,TG,CG,UG,XX,TT,DT,XA,XB,R1, 

R2.R3) 

This  subroutine  calculates  the  right  hand  side  of  the  gas 
continuity  (Rl),  momentum  (R2)  and  energy  (R3) , and  their 
derivatives  with  respect  to  p and  u. 

The  only  requirement  on  the  terms  in  the  call  list  other 
than  the  variables  used  in  Rl,  R2  and  R3 

K1  ■ 1 calculate  R1,R2,  and  R3 

2 " " ",  R3  and  their  derivatives. 

64.  SHKFRT(L,K1,K2,KK3,MAXIT, JJ1, JJ2,  JJ3) 

This  subroutine  controls  the  calculation  of  properties  in 
fron  of  and  behind  a shock  wave. 

The  terms  in  the  call  list  are: 
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K1 

K2 


KK3 

MAXIT 

JJ1 

JJ2 

JJ3 


number  of  the  shock 

0 normal  shock  front 

1 SHKFRT  called  from  BNSH  with  special  treatment  of 
A and  B properties 

2 not  used 

3 SHKFRT  called  from  BNSH  for  initial  fir9t  ,uess 
not  used 

maximum  number  of  iterations 
not  used 

ft  It 

ft  It 


65.  SHK1N (L , NR1L2 , N2 , M2 , W, X, Y , Z) 

This  subroutine  controls  the  initiation  of  a shock  wave. 

The  terms  in  the  call  list  are: 

NR1L2  - 1 right  traveling  shock  wave 
- 2 left  " " " 

N2  - the  point  number  of  the  point  behind  the 
inserted  shock  wave 

M2  ■ the  point  number  of  the  first  mesh  point 
in  front  of  the  shock. 

W ■ not  used 
X - " " 

y m M II 

Z - " " 

66.  SHKINT(L,K1, MAXIT, KK1,KK2,KK3,KK4) 

This  subroutine  controls  the  calculation  of  a shock  wave 
interacting  with  a gas-particle  interface. 

The  terms  in  the  call  list  are: 

K1  ■ number  of  the  shock  wave 
MAXIT  ■ maximum  number  of  iterations 


67.  SHOKEQ(L,MS,MQ) 

This  subroutine  solves  the  Rankin-Hugonlot  shock  relations 
for  the  properties  behind  a shock  wave  given  either  u,  p or  the 
shock  velocity,  U. 

The  terms  in  the  call  list  are: 

MS  ■ 1 right  traveling  shock 
2 left  " " 

MQ  - 1 specify  p benind  the  shock 

2 « u M ••  ii 

3 " the  shock  velocity,  II. 
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68.  SHTR(L,J,MM,I) 

This  subroutine  controls  the  calculation  of  properties  across 
a shock  (currently  not  in  use) 

The  terms  in  the  call  list  are: 

J - shock  number 
MM  ■ 1 right  traveling  shock 
2 left  " " 

I  ■ point  number  behind  the  shock 

69.  SMQ(L,  ISM, RW,EW,PW,PRW,PEW,NA,NB,XA,AB,TW) 

This  subroutine  calculates  the  pressure  PW,  temperature  TW 
and  the  derivatives  of  pressure  with  respect  to  density,  RW  and 
energy,  EW,  given  density  and  energy. 

The  terms  in  the  call  list  are: 

ISM  ■ 1 calculate  p 

2 calculate  p and  3p/3p 

3 calculate  p,  3p/3p  Celt!  op/  3E 

4 calculate  p and  3p/3E 

70.  SPPT(L,KKl,K2fKK3,PG,TG,RG,EG,m,XX2,XX3)  USER  SPECIFIED 

This  subroutine  calculates  C ,C  and  y 

v p 

The  terms  in  the  call  list  are: 

KK1  - not  used 

K2  - 1 Vlrlal  equation  of  state  is  used 

2 Van  der  Waal  equation  of  state  Is  used 

3 Ideal  gas  equation  of  state  us  used 

71.  SUPINF(Y,YAV,MAXY,MINY,K) 

This  subroutine  determines  the  maximum  MAXY , minimum  MI NY, 
and  average,  YAV  values  of  a variable  Y in  the  call  list,  and 
calculates  the  average  change  in  Y between  adjacent  mesh  points  on 
a given  time  line. 

The  terms  in  the  call  list  are: 

K - 1 behind  bullet 

2 in  front  of  bullet 
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72.  TBC0ND(L,K1,K2,K3,K4,XA,XB,UG) 

This  subroutine  specifies  the  boundary  conditions  for  the  gas 
are  present  at  the  boundary. 

call  list  are: 

boundary 
t " 
used 

ft 
«• 

• I 
fl 

1 

73.  WQ ( SG , SP , EG , EP , UG , UP , PG , PP , RG , RP , TG , TP , ES , XX , DT , TT , 

K1,K2,K3, ZZ, Z2,L) 

This  subroutine  calculates  the  rate  at  which  gas  mass  is 
created  per  unit  volume  of  mixture  during  burning,  w. 

The  only  requirement  on  the  terms  in  the  call  list  is  that 
all  variables  used  in  to  be  available. 

74.  WWGQ(SG,SP,EG,EP,UG,UP,PG,PP,RG,RP,TG,TP,ES,XX,DT,TT, 

K1,K2 ,K3, Zl, Z2, L) 

This  subroutine  calculates  the  rate  of  decrease  in  gas  mass 
per  unit  volume  of  mixture  due  to  losses  through  the  duct  wall, 

The  only  requirement  on  the  terms  in  the  call  list  is  that 

all  variables  used  in  u>  be  available. 

w 

75.  WWPQ  (SG , SP , EG , EP , UG , UP , PG , PP , RG , RP , TG , TP , ES , XX,  DT , TT , 
K1,K2,K3,Z1,Z2,L) 

This  subroutine  calculates  the  rate  of  decrease  in  particle 

mass  per  unit  volume  of  mixture  due  to  losses  through  the  duct 

wall,  to  . 
wp 

The  only  requirement  on  the  terms  in  the  call  list  is  that 

all  variables  used  in  to  be  available. 

wp 

76.  XINT(A,B,C,D,E) 

This  subroutine  performs  linear  interpolation  using  the 
formula 

XINT  - A + (B-A)E/D. 


phase  when  both  phases 
The  terms  in  the 


K1 

K2 

K3 

K4 

XA 

XB 


1 left 

2 righ 
not 
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77.  VSWICH(L,K1,K2,K3,K4,K5,K6,KK7,KK8,K9,K10) 

This  subroutine  controls  the  switching  pf  properties 
from  one  location  to  another. 

The  terms  in  the  call  list  are: 

K1  - first  point  to  be  switched 

K2  - last  " ” " " 

K3  ■ 1 always 

K4,K5  ■ properties  will  be  switched  from  line  K5 
to  line  K4. 

KG  - 1 K7  and  K8  ■ I (K1  < I < K2) 

- 2 K7  and  K8  - KK7  and  KX8 
K7,K8  - point  K8  will  be  switched  to  K7  location 
KK7.KK8  - (see  K6  - 2) 

K9  - 1 always 
K10  - 0 always 

III.  Description  of  the  options  available  for  a user  which  re- 
quire subroutines  to  be  altered 

1.  Change  the  equation  of  state 

a.  To  change  the  equation  of  state  the  user  must  put  the 
equation  of  state  in  the  following  form: 

p • P(P.T) 

and  he  must  develope  relations  of  the  form: 
f(E,p,T)  - 0 

and 

g(S,p, T)  - 0 

b.  Subroutines  PTQ,  QET,  NTROPY  and  FRET  must  be  changed 
in  accordance  with  the  Instructions  listed  In  the 
program. 

2.  Change  the  expression  for  the  regression  rate 

The  user  must  alter  subroutine  DZDTQ  In  accordance 
with  instructions  In  the  program. 
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J.  Change  the  shape  of  the  particles. 

The  user  must  change  the  expressions  for  the  lolume  of  the 
particle  VP  as  a function  of  regression  distance  zz  and  its 
time  derivative  in  subroutine  QGQ  and  WGQ. 

4.  Change  the  drag  force  expression. 

The  user  must  alter  subroutine  FQ  in  accordance  with 
instructions  listed  in  the  program. 

5.  Change  the  expression  for  heat  transfer  with  the  wall. 

The  user  must  alter  subroutine  QWGQ  for  the  gas  phase  and 
QWPQ  for  the  particle  phase  in  accordance  with  instructions 
listed  in  the  program. 

6.  Change  the  expression  for  mass  transfer  with  the  wall. 

The  user  must  alter  subroutine  WWGQ  for  the  gas  ohast  and 
WWPQ  for  the  particle  phase  in  accordance  with  instructions 
listed  in  the  program. 

7.  Change  the  expression  for  the  friction  force  at  the  wall. 

The  user  must  alter  subroutine  FWGQ  for  the  gas  phase  and 
FWPQ  for  the  particle  phase  in  accordance  with  instructions 
listed  in  the  program. 

8.  Change  the  expression  for  the  area  of  the  duct  as  a 
function  of  x. 

The  user  must  alter  subroutines  AR  and  ARX  in  accordance 
with  the  instructions  listed  in  the  program. 

9.  Change  the  expression  for  C^. 

The  user  must  alter  subroutine  SPPT  by  changing  the  expression 
for  CV. 
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IV  Instructions  for  running  the  computer  code  TWOFLO  (additional 

information  can  he  found  in  the  sample  output,  Sec. (6)) 

1.  Before  proceeding  to  discuss  the  input  cards  necessary  to  run 
TWOFLO,  let  us  first  explain  the  restart  option  available  to  the 
user.  This  program  is  designed  to  allow  the  user  to  divide  a long 
run  up  into  several  shorter  ones,  or  to  extend  a run  that  has  been 
terminated  due  to  process  time  limits.  The  important  point  to 
he  noted  is  that  the  data  is  dumped  out  by  subroutine  AOUT  and 
read  back  into  the  program  in  Identical  format  in  subroutine  AIN. 

If  one  wishes  to  change  to  a form  of  output  other  than  a dump  file, 
(ic,  punch  cards,  etc.)  his  only  concern  need  be  that  all  data 
currently  being  transferred  out  of  the  program  in  the  present  form 
be  transfered  out  in  the  new  form.  This  data  must  then  be  read 
back  into  the  program  in  entirety. 

The  procedure  for  utilizing  the  restart  capability  when 
running  TWOFLO  is  as  follows: 

a.  Define  three  disk  or  tape  files. 

b.  Set  KDUMPN  ■ 1 (data  is  read  in  from  card  deck)  for  initial 
run  and  KDUMPN  * 2 (data  is  read  in  from  a file)  for  subsequent 
restarts. 

c.  Set  KDUMPT  * 2 (data  will  be  dumped  on  file  for  all  time 
lines  evenly  divisible  by  NOUTF) . 

d.  Set  the  file  codes  IFILEI,  IFILE1,  IFILE2  and  IFILE3. 

The  data  will  be  read  into  the  program  from  IFILEI  and  immediately 
dumped  on  IFILE3  for  storage.  Subsequent  time  lines  are  dumped 
on  both  IFILE2  and  IFILE3  in  accordance  with  NOUTF. 


2.  The  following  is  a detailed  description  of  the  data  cards 
needed  to  initiate  a calculation  using  TWOFLO. 


VARIABLE 

NAME 


KTLOT 


N PLOTS 


NPLOTP 


FORMAT 


9-12  14 


0 (no  plots  will  be  made) 

1 (plots  of  prescribed 
variables  vs.  position 
will  be  made. 

no.  of  variables  to  be 
plotted.  It  is  possible 
to  plot  up  to  six  variables 
vs.  position. 

Every  time  line  divisible 
by  NPLOTP  will  be  plotted. 


If  KPLOT  - 0 omit  cards  no.  2,3  and  4 


IPL0T6 


NOADDB  5 1-4 


NOADDA 


■ 0 (No  plotting) 

plots  will  be  made 


Code  indicating  which  var- 
iables are  to  be  plotted. 
(See  subroutine  INIT  for 
variable  codes) 


Nos.  of  plots  of  a partic- 
ular variable  per  graph 


0 Program  will  automatically 
add  points  behind  the 
bullet  if  required. 

1 No  pts.  will  be  added  be- 
hind the  b'  llet 

0 Program  will  automatically 
add  points  in  front  of  the 
bullet  if  required. 

1 No  pts.  will  be  added  in 
front  of  the  bullet. 

DTMAX  before  interface  is 
NOADDI  times  as  large  as 
normal 


Maximum  no.  of  pts.  allowed 
behind  the  bullet 

Maximum  no.  of  pts.  allowed 
in  front  of  bullet. 


VARIABLE 

NAME 


CARD 

No. 


COLUMN 

No. 


FORMAT 


DESCRIPTION 


PTES 

7 

1-15 

E15.8 

Limiting  value  for  volume 
fraction  of  particles,  e.  If 
c at  any  point  is  less  than 
PTES  particle  effects  are 
Ignored. 

NTCAL 

8 

1-4 

14 

No.  of  special  times  lines 
at  which  the  properties  will 
be  Interpolated  for  and 
printed  out.  (up  to  five 
special  times  allowable) . 

If  NTCAL  - 0 omit  card 

No.  9 

TCAL(I) 

9 

1-75 

5E15.8 

The  special  times  at  which 

I- 1, NTCAL 

the  properties  will  be 
printed  out. 

KDUMPN 

10 

1 

11 

- 1 

(If  KDUMPN  - 2 data  is  to  be 
read  in  from  existing  disk 
or  tape  file.  See  next 
section  for  explanation  of 
input) 

KDUMPT 

2 

11 

- 1 

if  calculations  are  not  to 
be  dumped  on  a file 

- 2 

if  calculations  are  to  be 
dumped  on  existing  disk  or 
tape  file. 

IFILE1 

3 

11 

If  KDUMPT  - 1,  the  values 

IFILE2 

A 

11 

of  IFILEI  and  IFILE2  are 

meaningless  (set  equal  to 
zero) 

If  KDUMPT  - 2 calculated 
time  lines  are  dumped  on 
disk  or  tape  files  nos. 
1FILE2  and  IFILE2. 

IFILE3 

5 

11 

If  KDUMPN  - 1 value  of 
IFILE3  is  meaningless  (set 
equal  to  zero) . 

If  KDUMPN  ■ 2 the  calcula- 
tions read  in  from  IFILEI 
are  dumped  on  IFILE3  as  an 
auxiliary  file 

IFILEI 

5 

11 

If  KDUMPN  - 2 Initial  time 
line  is  read  in  from  IFILEI. 
If  KDUMPN  - 1 value  of  IFILEI 

is  meaningless  (set  equal  to 
zero) 

12 


NOUTF 


NOUTP 


I DRAG 


Any  time  line  number  into 
which  NOUTF  divides  evenly 
will  be  dumped  on  IFILE1  and 
1FILE2. 

Any  time  line  into  which 
NOUTP  divides  evenly  will  be 
rinted  out. 


B-CONSTANTSi  13-42  1-60  4E15.8 


PERPGD 

PERUGD 

PERRGD 


PERPGA 

PERUGA 

PERRGA 


ITMAX 


31-45  E15.8 


44  1-15 


E15.8 


45  1-15  E15.8 


16-19 


20-23 


24-38  E15.8 


On  the  next  thirty (30)  cards 
the  B-CONSTANTS  are  read  in, 
four(4)  values  per  card  (a 
total  of  120  CONSTANTS) 

Card  No.  eleven (11)  would 
contain  B(l,l),  B(l,2), 
B(l,3),  and  B(l,4)  card 
twelve  B(l,5),  B(l,6),  B(l,7) 
B(l, 8)  card  No.  26  B(2,l), 
B(2, 2)  B(2, 3) , B(2,4)  etc. 

(a  description  of  the  b 
CONSTANTS  appears  at  the  end 
of  this  section)  . 

. 50000000E-02 

. 50000000 E-02 

. 50000000E-02 

. 10000000E-02 
. 10000000E-02 

. 10000000E-02 

The  time  of  the  initial  time 
line 

Time  line  number  of  the 
initial  time  line 

No.  of  time  lines  to  be 
calculated 

Time  of  last  time  line  cal- 
culated will  be  less  than 
TMAX  (This  parameter  over- 
rides NDELT) 
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DTMIN 


46 


1-15 


E15.8 


1-4 

14 

5-8 

14 

9-12 

14 

13-16 

14 

17-20 

14 

1-15 

E15.8 

16-30 

E15.8 

31-45 

E15.8 

46-60 

E15.8 

61-75 

E15.8 

1-15 

16-30 

31-45 

46-60 

61-75 


E15.8 
E15.8 
Elc.  8 


The  minimum  value  for  DELTA  t 
(time) 

The  maximum  value  for  DELTA  t 


If  calculated  value  of  DT  is 
greater  than  DTFIX,  DT  will 
be  set  equal  to  DTFIX  without 
adding  points 

Maximum  number  of  iterations 

If  pressure  gradient  between 
two  successive  points  on  time 
line  is  greater  than  PXTOL 
and  ISMAX  > 0,  a shock  1* 
Inserted 


0 number  of  area  discontinui- 
ties (not  including  gas- 
particle  interface)  in  the 
problem 


The  point  numbers  of  the 
area  discontinuities 

Point  number  of  breech 

Point  number  of  left  surface 
of  bullet 

Point  number  of  Interface 
between  gas  and  two-phase 
region 

Point  number  of  right  surface 
of  bullet 

Point  number  of  muzzle 

The  following  twelve  variable! 
are  the  tolerances  used  in  th< 
program's  iterating  procedure 
e.g.,  TOLSG  is  the  tolerance 
for  SG(gas  concentration)  ee 
output  for  further  informatioi 


TOL.ES 

54 

1-15 

E15.8 

TOXX 

16-30 

E15.8 

MEOS ( 1 ) 

55 

1-2 

12 

MEOS (2) 


ISMAX 


If  I! 

IS  (1> 
ISRIL2(1) 


56 

1-4 

14 

MAX  - 0 

omit  cards  No.  57 

57 

1-4 

14 

5-8 

14 

9-23 

E15.8 

PGSH(l)  | | 9-23  I E15.8 

If  ISMAX  ■ 1 omit  card  No.  58 


* 1 A virial  equation  of  state 
is  used  in  region  behind  the 
bullet 

■ 2 Van  der  Waals  equation  of 

state  used  in  region  behind 
the  bullet 

■ 3 Noble-Abel  equation  of  state 

is  used  in  region  behind  the 
bullet 

■ 1,2  or  3 same  as  ME0S(1) 

except  region  is  in  front  of 
the  bullet. 

» number  of  shocks  initially 
Inserted  (<2) 


point  number  where  shock 
number  1 is  to  be  inserted. 

1 right  traveling  shock  is 
inserted 

2 left  traveling  shock  is 
inserted 

pressure  behind  the  shock 


IS(2) 


ISR1L2(2) 


PGSH(2) 

LCODE 


1-4 


9-23  E15.8 


59  1-4 


point  number  where  shock 
number  2 is  to  be  inserted 
(must  be  in  front  of  bullet) 

1 right  traveling  shock  *s 
inserted. 

2 left  traveling  shock  is 
inserted 

pressure  behind  the  shock 

1 Properties  of  all  pts.  behin< 
bullet  are  initially  the 
same.  There  will  be  IFRE(2) 
-1  such  pts.,  each  separated 
by  a DELTA  X of  DX  (calcu- 
lated in  subroutine  INIT) 
only  one  set  of  properties 
behind  bullet  should  be  read 
in 
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LCODEB 


0 Properties  of  all  pts.  be- 
hind bullet  must  be  read  in 

1 Properties  of  all  pts.  in 
front  of  the  bullet  are 
initially  the  same.  There 
will  be  IFRE(4)-IFRE(3)-1 
such  pts.,  each  being 
separated  by  a DELTA  X of 
(calculated  in  INIT)  only 
one  set  of  properties  in 
front  of  bullet  should  be 
read  in. 

0 Properties  of  all  pts.  in 
front  of  bullet  must  be 
read  in. 


If  LCODE  - 1 

i - 2 

L 

one  set 
IFRE(2) 

of  cards 
sets  of  c 

IK 

60 

1-4 

14 

IJ 

5-8 

14 

XX(l.I) 

9-23 

E15.8 

UG(l.I) 

24-38 

E15.8 

RG(l.I) 

39-53 

PG(l.I) 

54-68 

E15.8 

Point  number  of  properties 
on  initial  time  line 

Distance  from  IFRE(l)  on 
initial  time  line. 


Gas  velocity 
Gas  density 
Gas  pressure 


IJ  > INT(l)  do  not  read  cards  61  and  62 


E1S.8 

E15.8 

E15.8 


ZZ(l.I) 


Particle  velocity 

0.0 

Particle  density 

Volume  fraction  of  particles 


Regression  distance 


If  LCODEB  - 1 and  IFRE(4)  + IFRE(2)  read  card  63 

- 2 and  " " " " " " ” " IFRE(4)- 

IFRE(2)  time 


Same  as 

card  60 

63 

■ 1 Causes  nondimensionalization 

■ 0 No  nondimensionalization 


3.  The  following  data  cards  are  used  to  restart  TWOFLO 


CuRD 

COLUMN 

FORMAT  1 

DESCRIPTION 

No. 

No. 

The  f 

ollowing  data  cards  are  used  to  restart  the  program. 

Cards 

1-9 

are  the  same  as  the  Initial  data 

Card 

10 

is  the  same  as  the 

Initial  data  except  KDUMPN-2 

Cards 

11-12 

are  the  seme  as  the  initial  data 

Cards 

13-15 

are  the  same  as  card  52-54  in  the  initial  data 

T!1AX 

16 

1-15 

E15.8 

time  of  last  time  line  cal- 

culated  will  be  less  than 

TMAX  (this  parameter  over- 

rides  NDELT) . 

E15.8 


18-32  E15.8 


33-36 


37-40 


1 Value  of  DELTA  T remains  the 
same  (no  change  for  restart) 

2 Value  of  DELTA  T is  set 
equal  to  DTNEW 

If  NDTNEW  - 1,  value  of  DTNE 
is  meaningless  (set  equal  to 
zero) 

If  NDTNEW  - 2,  DELTA  T is 
set  equal  to  DTNEW 

1 Value  of  DTMIN  (see  previous 
section  for  explanation) 
remains  same 

2 DTMIN  is  set  equal  to  DTMINW 

If  NDTMIN  - 1,  value  of 
DTMINW  is  meaningless  (set 
equal  to  zero) 

If  NDTMIN  - 2,  DTMIN  is  set 
equal  to  DTMINW 
The  number  of  time  lines  to 
be  calculated  on  this  run. 

1 Value  of  DTFIX  remains  same. 

2 Value  of  DTFIX  is  set  equal 
to  DTFIXW 


DTFIXW 


NDTMAX 


DTMAXW 


NINT 


INTW 


NPXTOL 


PXTOLW 


41-55  E15.8  If  NDTFIX  « 1,  value  of 

DTFIXW  is  meaningless  (set 
equal  to  zero) 

If  NDTFIX  - 2,  DTFIX  is  set 
equal  to  DTFIXW. 

18  1 II  ■ 1,  Value  of  DTMAX  remains  the 

same. 

■ 2 Value  of  DTMAX  is  changed 

to  DTMAXW 

2-16  E15.8  If  NDTMAX  - 1,  value  of 

DTMAXW  is  meaningless  (set 
equal  to  zero) 

If  NDTMAX  - 2,  DTMAX  is  set 
equal  to  DTMAXW 

19  1 II  ■ 1,  Value  of  INT(l)  remains  the 

same 

■ 2,  Value  of  INT(2)  is  set  equal 

to  INTW 

2-4  13  If  NINT  - 1,  Value  of  INTW 

is  meaningless  (set  equal  to 
zero) 

If  NINT  - 2,  Value  of  INT(l) 
is  set  equal  to  INTW 

20  1-4  14  - 1,  Value  of  PXTOL  remains  the 

same 

* 2,  Value  of  PXTOL  is  set  equal 
to  PXTOLW  (If  e is  less 
than  PXTOL  particle  proper- 
ties are  ignored) 

5-19  E15.8  If  NPXTOL  - 1 value  of 

PXTOLW  is  meaningless  (set 
equal  to  zero) 

If  NPXTOL  - 2,  PXTOL  is  set 
equal  to  PXTOLW. 


T39PW 


4.  Description  of  the  B(2,60)  array  input  constants. 


The  "B"  constants  which  specify  material  properties  will  be 
underlined.  Such  properties  designated  by  the  1 location  in  the 
array  (ie.  B(l,  ))  are  for  the  material  behind  the  bullet;  while 
those  designated  by  the  2 location  are  for  the  material  in  front 
of  the  bullet.  Thus  it  is  only  necessary  to  explain  the  B(l,  ) 
portion  of  the  array.  Any  constants  that  are  not  explained  below 
are  currently  not  defined  in  the  code  and  can  be  set  equal  to 
zero. 


B(M) 

B(1 ,2) 
B(l,3) 
B(l,4) 
B( 1 , 5) 

B( 1 ,6) 

B ( 1 ,8) 
B(l.ll) 
B(1 , 12) 
B(l, 13) 
B(l,14) 
B(1 , 17) 
B(1 ,20) 
B ( 1 ,21) 
B(l,22) 
B ( 1 ,23) 
B(1 ,24) 
B ( 1 ,25) 

B ( 1 ,26) 
B(1 ,27) 
B ( 1 ,28) 
B ( 2 , 29) 
B(1 , 30) 
B(1 , 31) 
B ( 1 , 32) 
B ( 1 , 33) 
BC1.34) 
B(l,35) 
B( 1 , 36) 

ILUIZI 
B(1 »38) 
B(l,40) 
B(l,41) 
B ( 1 , 42 ) 
B(i ,43) 


B , 


Reference  density  for  Equation  of  State. 

Reference  temperature  for  Equation  of  State. 
Reference  energy  for  the  Equation  of  State. 
Reference  pressure  for  Equation  of  State. 
Geometrical  parameter  used  in  subroutines  QGQ 
and  WGQ  to  describe  the  particle  shape. 
Geometrical  parameter  used  in  subroutines  QGQ 
and  WGQ  to  describe  the  particle  shape. 

Hold  back  pressure 
Universal  gas  constant 
Mass  of  gas  per  mole 

Constant  in  Van  der  Waal  Equation  of  State 

Constant  in  Van  der  Waal  Equation  of  State 

Wall  temperature 

Nondimens ionalizing  density 

Nondimensionalizing  distance 

Nondimens ionalizing  pressure 

Initial  x location  of  the  bullet 

Bullet  length 

, Barrel  length  (measured  from  the  beginning  of  the 
barrel  diameter) 

Diameter  of  the  chamber 
x location  where  chamber  diameter  ends 
Barrel  diameter 

x location  of  the  beginning  of  the  barrel  diameter 
Reference  dynamic  viscosity  used  to  calculate  u 
Reference  temperature  used  to  calculate  u 
Reference  pressure  used  to  calculate  u 
Constant  used  to  calculate  u 
Constant  used  to  calculate  p 

Particle  velocity  ratio  across  area  discontinuity 
Reference  temperature  used  in  entropy  equation 
Constant  in  equation  for  C 
Constant  In  equation  for  Cv 

Specific  internal  energy  o?  combustion  products 
Heat  released  per  unit  volume  of  propellant  burnt 
Constant  in  regression  formula 
Constant  in  regression  formula 
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B( 1 , 48) 
B ( 1 ,49) 
B ( 1 , 50) 
B ( 1 ,51) 
B(1 , 52) 

B(1 ,53) 


F^,  Bore  resistance  force 
Ag,  Cross-sectional  area  of 
Mg,  Bullet  mass 
PM,  Muzzle  pressure 
tgD>  The  time  that  it  takes 
after  the  shock  leaves 
tfip.  The  time  that  it  takes 
after  the  bullet  leaves 


the  bullet 


to  reach  p 

JKi  M 


for  p 
the  b 

for  p to  reach  p 
the  barrel  M 


5.  Description  of  user  specified  functions  and  variables  currently  in  TWOFLO. 

a.  Drag  Force  (spherical  particle) 


F “ - CD  • e • p • lu-uj  • (u-uJ/B(l,5) 


where 

CD  - ||  + 4-RE(_l/3) 
RE 


RE  - 2*p*B(l,5)*|u-u  |/u 

U • B(l,30)  + [t-B(1,31)] -B(l,33)  + [p-B(l ,32)]  *B(1 ,34) 

b.  Constant  volume  specific  heat  Cv 

Cv»  B(L,37)  + B(L,38)-T  (L-1,2) 

c.  Equation  of  State  (Ideal  gas) 

p - B(L, 11) *p *T/B(L , 12)  (L-1,2) 

d.  Volume  of  a particle  (pancake) 

V - ir*  [B (1,5)  — Z j 2*{B(l,6)+4*  [b(1,5)-Z]  /3} 

e.  various  parameters  in  Eq.  (3) 

4>  • Fu  /o 
P 

u - u 
w 

v - 0 

w 

E - E 
w 
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h.  Sample  initial  data 

The  sample  input  data  presenced  here  correspondes  to  computer 
runs  A,B  and  C presented  in  the  Results  and  Conclusions  section  of 
the  main  report. 

CARD  1 KPLOT  « 0 

NPLOTS  * 0 

NPLOTP  - 0 

* 

CARD  2 Not  used  (omit  this  card) 

CARD  3 Not  used  (omit  this  card) 

CARD  4 Not  used  (omit  this  card) 

CARD  5 NOADDB  = 0 

NOADDA  - 1 

NO ADD I = 3 

CARD  6 MAXPT  - 30 

MAXPTB  * Not  used  (an  arbitrary  value  of  0 is  read  in) 

CARD  7 PTES  - 0. 10000000E-02 

CARD  8 NTCAL  - 4 

CARD  9 TCAL(l)  - 0. 20000000E-03  s 

TCAL(2)  - 0. 40000000E-03  s 

TCAL(3)  - 0. 60000000E-03  s 

TCAL(4)  - 0 . 80000000E-03  s 

CARD  10  KDUMPN  * 1 

KDUMPT  - 2 

IFILE1  - 1 

IFILE2  - 2 

I FILE 3 - 1 

IFILEI  - 1 

CARD  11  NOUTF  - 1 

NOUTP  - 1 

CARD  12  IDRAG  - 1 

CARD  13  B( 1 , 1)  - 0. 12800000E+01  kg/m3 

B(  1 , 2)  ■ Not  used  (an  arbitrary  value  of  0.0  K is  read  in) 

B(l,3)  ■ Not  used  (an  arbitrary  value  of  0.0  J/kg  is  reed  in) 

B(l,4)  ■ Not  used  (an  arbitrary  value  of  0.0  Pa  is  read  in) 

CARD  14  B(l,5)  • 0.27305000E-03  o (Runs  A & C)  0. 13650000E-03  m (Run  C) 

8(1,6)  - 0. 38100000E-03  m (Runs  A & C)  0. 19050000E-03  m (Run  C) 

B(l,7)  ■ Not  defined** 

B(l,8)  - 0.20700000E408  Pa 

* The  term  Not  used  indicates  that  the  variable  is  defined  but  not  used 
for  this  specific  run. 

**The  term  Not  defined  Indicates  that  the  term  is  never  used  in  the  program. 

Ill 


t 


CARO 


CARD 


CARD 


CARD 


CARD 


CARD 


CARD 


CARD 


CARD 


CARD 


15  B( 1 ,9) 
B(l,10) 
B(l,ll) 
B(1 , 12) 


Not  defined  = 0 . OOOOOOOOE+OO 
Not  defined  * 0 . OOOOOOOOE+OO 
0. 83140000E+01  J/(K-rool) 

0. 25547800E-01  kg/tnol 


16  B( 1 , 13) 

B( 1 , 14) 
B(1 , 15) 
B( 1 , 16) 


= Not  used  (an  arbitrary  value  of  0.0  N*m^/kg2  is 
read  in) 

*■  Not  used  (an  arbitrary  value  of  0.0  mVkg  is  read  in) 
= Not  defined  ■ 0 . OOOOOOOOE+OO 
* Not  defined  ■ 0 . OOOOOOOOE+OO 


17  B(1 , 17) 
B(1 , 18) 
B(1 , 19) 
B(1 ,20) 


0. OOOOOOOOE+OO  K 
Not  defined  ■ 0. OOOOOOOOE+OO 
Not  defined  ■ 0. OOOOOOOOE+OO 
0. 15000000E+03  kg/m3 


18  B(1 ,21) 

B ( 1 , 22) 
B(1 ,23) 
B ( 1 ,24) 


0.40000000E-01  m 
0. 20000000E+09  Pa 
0. 33100000E-01  m 
0. 15000000E-01  m 


19  B ( 1 , 25 ) 
B(1 ,26) 
B(1 ,27) 
B(1 ,28) 


0. 47000000E+00  tn 
0.90750000E-02  m 
0.26850000E-01  m 
0. 56500000E-02  m 


20  B(1 ,29) 
B(1 ,30) 
B ( 1 , 31 ) 
B(1 ,32) 


0.31  OOOOOOE-O 1 tn 
0. 18192000E-04  N-a/tn2 
0.29300000E+03  K 
0. 10130000E+06  Pa 


21  B ( 1 ,33) 
B(l,34) 
B(1 ,35) 
B(l,36) 


0. 53600000E-07  N-s/(m2.K) 

0. 12080000E-12  s 

Not  defined  ■ 0 . OOOOOOOOE+OO 

0. 29800000E+03  K 


22  B ( 1 ,37) 
B( 1 , 38) 
B ( 1 , 39) 
B(1 ,40) 


0. 13559000E+04  J/(kg-K) 

0. OOOOOOOOE+OO  J/(kg-K2) 

Not  defined  * 0 . OOOOOOOOE+OO 
0. OOOOOOOOE+OO  J/kg 


23  B ( 1 , 4 1 ) 
B(1 ,42) 
B ( 1 ,43) 
B(l,44) 


0. 63000000E+10  J/m3 

0. 12000000E-07  m/(Pa6-a)  (6-B(l,43)) 

0. 84550000E+00 

Not  defined  ■ 0. OOOOOOOOE+OO 


24  B ( 1 ,45) 
B(1 ,46) 
B ( 1 ,47) 
B(1 ,48) 


Not  defined  ■ 0 . OOOOOOOOE+OO 
Not  defined  ■ 0 . OOOOOOOOE+OO 
Not  defined  ■ 0 . OOOOOOOOE+OO 
0. 16700000E+03  N 
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CARD  25  R( 1 ,49)  = 0 . 24279500E-04  m2 

B ( 1 , 50)  * 0. 36290000E-02  kg 

B ( 1 ,51)  = 0. 10132500E+06  Pa 

B ( 1 ,52)  = 0. 20000000E-04  s 

CARD  26  B(1 , 53)  = 0. 50000000E-04  s 

B(  1 ,54)  = Not  defined  - 0 . 00000000E+00 

B(l,5b)  * Not  defined  - 0 . OOOOOOOOE+OO 

B ( 1 , 56)  * Not  defined  - 0 . 00000000E+00 

CARD  27  B(  1 ,57)  - B(l,60)  - Not  defined  - 0. OOOOOOOOE+OO 

CARD  28  B(2,l)  * Not  used  (an  arbitrary  value  of  0.0  kg/m^  is  read  in) 

B(2,2)  * Not  used  (an  arbirtary  value  of  0.0  K is  read  in) 

B(2,3)  * Not  used  (an  arbitrary  value  of  0.0  J/kg  is  read  in) 

B(2,4)  * Not  used  (an  arbitrary  value  of  0.0  Pa  isread  in) 

CARD  29  B(2 , 5)  - B(2,8)  - Not  defined  - 0. OOOOOOOOE+OO 

CARD  30  B(2,9)  - Not  defined  - 0 . 00000000E+00 

B(2, 10)  * Not  defined-  0 . OOOOOOOOE+OO 

B(2,ll)  - Not  used  (an  arbitrary  value  of  0.0  J/(K*mol)  is  read  in) 

B(2,12)  - Not  used  (an  arbitrary  value  of  0.0  kg/mol  is  read  in) 

CARD  31  B (2 , 13)  - Not  used  (an  arbirtary  value  of  0.0  N*m6/kg2  is 

read  in) 

B(2,14)  - Not  used  (an  arbitrary  value  of  0.0  m^/kg  isread  in) 

B(2, 15)  - Not  defined  - 0 . OOOOOOOOE+OO 

B(2 , 16)  - Not  defined  - 0 . OOOOOOOOE+OO 

CARD  32  B(2, 17)  - B(2,20)  - Not  defined  - 0. OOOOOOOOE+OO 

CARD  33  B(2 ,21)  - B(2,24)  - Not  defined  - 0 . OOOOOOOOE+OO 

CARD  34  B(2,25)  - B(2,28)  - Not  defined  - 0. OOOOOOOOE+OO 

CARD  35  B(2,29)  - B(2,32)  - Not  defined  - 0. OOOOOOOOE+OO 

CARD  36  B(2,33)  - Not  defined  - 0 . OOOOOOOOE+OO 

B(2,34)  • Not  defined  - 0 . OOOOOOOOE+OO 

B(2,35)  - Not  defined  - 0 . OOOOOOOOE+OO 

B(2,36)  - Not  used  (an  arbirtary  value  of  0.0  K is  read  in) 

CARD  37  B(2,37)  - Not  used  (an  arbitrary  value  of  0.0  J/(kg*K)  is  read  in) 

B(2,38)  - Not  used  (an  arbitrary  value  of  0.0  J/(kg*K2)  is  read  in) 

B(2, 39)  - Not  defined  - 0 . OOOOOOOOE+OO 

B(2,40)  - Not  defined  - 0 . OOOOOOOOE+OO 
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CARD  38 
CARD  39 
CARD  40 
CARD  41 
CARD  42 
CARD  43 

CARD  44 


B(2,41) 

B(2 ,45) 

B(2,49) 

B(2,53) 

B(2 , 57) 

PERPCD 

PERUGD 

PERRGD 


B(2,44) 
B(2 , 48) 
B(2 , 52) 
B(2 , 56) 
B(2 ,60) 


Not  defined 
Not  defined 
Not  defined 
Not  defined 
Not  defined 


O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO 


0.50000000E-02 
0. 50000000E-02 
0. 50000000E-02 


PERPGA 

PERUGA 

PERRGA 


0. 10000000E-01 
0. 10000000E-01 
0. 10000000E-01 


CARD  45  TIME 
NTIME 
NDELT 
TMAX 


O.OOOOOOOOE+OO  a 
0 
40 

0. 12000000E-02  s 


CARD  46 

CARD  47 
CARD  48 
CARD  49 

CARD  50 


DTMIN 

DTMAX 

DTFIX 

MAXIT 

PXTOL 

NDISC 

Not  used 


- 0. 45000000E-05  s 

* Not  used  (a  value  less  than  2.2  x DTMIN  should 
be  read  in)  ■ O.OOOOOOOOE+OO  s 

■ Not  used  ( a value  much  larger  than  DTMIN  should 
be  read  in)  « 0. 10000000E+10  s 

- 15 

- 0. 50000000E+01  Pa/m 

* 0 (Program  is  currently  not  operational  for 

NDISC  greater  than  0) 

(omit  this  card) 


CARD  51  IFRE(l)  - 1 

IFRE(2)  -=  6 

INT(l)  - 6 

IFRE(3)  - 6 

IFRE(4)  - 6 


CARD  52  TOLSG 
TOLEG 
TOLUG 
TOLPG 
TOLRG 


0. 10000000E-02 
0.10000000E-02 
0. 10000000E-02 
0. 10000000E-02 
0. 10000000E-02 
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s 


CARP 

S3 

TOLSP 

m 

0. 10000000E-02 

TOLEP 

St 

Not  used  (an  arbitrary  value 
Is  read  in) 

of 

0. 10000000E-02 

TOLUP 

- 

0. 10000000E-02 

TOLPP 

* 

Not  used  (an  arbitrary  value 
Is  read  in) 

of 

0. 10000000E-02 

TOLRP 

’ 

Not  used  (an  arbitrary  value 
is  read  in) 

of 

0. 10000000E-02 

CARD 

54 

TOLES 

= 

0. 10000000E-02 

TOLXX 

= 

0. 10000000E-02 

CARD 

55 

ME0S(1) 

m 

3 

MEOS(2) 

= 

Not  used  (an  arbitrary  value 

of 

3 is  read  in) 

CARD 

56 

ISMAX 

= 

0 

CARD 

57 

Not  used 

(omit  this  card) 

CARD 

58 

Not  used 

(omit  this  card) 

CARD 

59 

LCODE 

- 

1 

L.C0DEB 

= 

Not  used  (an  arbitrary  value 

of 

1 is  read  in) 

CARD 

60 

IK 

m 

1 

IJ 

X 

1 

XX(1,1) 

=■ 

O.OOOOOOOOE+OO  m 

UG(1,1) 

m 

O.OOOOOOOOE+OO  m/s 

RG(I,1) 

m 

m 

0. 10755000E+00  kg/m3  (Runs  A 
0.14637200E4-02  kg/m3  (Run  C) 

& 

B) 

PG(1 , 1) 

m 

0. 10132500E+06  Pa  (Runs  A & B) 

* 

0.13790000E+08  Pa  (Run  C) 

CARD 

61 

UP(1,1) 

m 

0 . 00000000E-H10  m/s 

EP(1,I) 

m 

Not  used  (an  arbitrary  value 
is  read  in) 

of 

O.OOOOOOOOE+OO  J/kg 

RP ( 1 , 1 ) 

m 

0. 16050000E-K14  kg/m3 

ES(1,1) 

m 

0. 57500000E+00 

CARD 

62 

ZZ(I,I) 

- 

O.OOOOOOOOE+OO  m 

CARD 

63 

Not  used 

(omit  this  card) 

CARD 

64 

NODIM 

m 

0 
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r«  SV*'  ‘ 


TOLPP 

* 

Not  used  (an  arbitrary  value 
is  read  in) 

of 

0. 10000000E-02 

TOLRP 

Not  used  (an  arbitrary  value 
is  read  in) 

of 

0. 10000000E-02 

CARD 

15 

TOLES 

B 

0. 10000000E-02 

TOLXX 

s 

0. 10000000E-02 

CARD 

16 

TMAX 

m 

0. 12000000E-02  s 

CARD 

17 

NDTNEW 

B 

1 

DTNEW 

B 

Not  used  (an  arbitrary  value 
is  read  in) 

of 

O.OOOOOOOOE+OO 

s 

NDTMTN 

B 

1 

DTMINW 

* 

Not  used  (an  arbitrary  value 
is  read  in) 

of 

0.00000000E+00 

5 

NDELT 

B 

5 

NDTFIX 

B 

1 

DTFIXW 

m 

Not  used  (an  arbitrary  value 
is  read  in) 

of 

O.OOOOOOOOE+OO 

S 

CARD 

18 

NDTMAX 

m 

1 

DTMAXW 

m 

Not  used  (an  arbitrary  value 
is  read  in) 

of 

0.00000000E+00 

s 

CARD 

19 

NINT 

B 

1 

INTW 

B 

Not  used  (an  arbitrary  value 

of 

0 is  read  in) 

CARD 

20 

NPXTOL 

m 

1 

PXTOLW 

m 

Not  used  (an  arbitrary  value  of 

O.OOOOOOOOE+OO 

Pa 

is  read  in) 
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ft  . 


7.  Sample  restart  data 


CARD 

1 

KPI.0T 

m 

0 

: 

N PLOTS 

- 

0 

NPLOTP 

- 

0 

CARD 

2 

Not  used 

(omit  this  card) 

CARD 

3 

Not  used 

(omit  this  card) 

CARD 

4 

Not  used 

(omit  this  card) 

CARD 

5 

NOADDB 

m 

0 

1 

I 

NOADDA 

B 

1 

i 

NOADDI 

B 

3 

CARD 

6 

MAXPT 

X 

30 

MAXPTB 

m 

Not  used  (an  arbitrary  value  of 

0 is  read  in) 

CARD 

7 

PTES 

m 

0. 10000000E-02 

I 

CARD 

8 

NT  CAL 

m 

4 

CARD 

9 

TCAL(l) 

m 

0. 20000000E-03  s 

TCAL(2) 

m 

0.40000000E-03  s 

TCAL(3) 

m 

O.60000000E-03  a 

TCAL(4) 

m 

0. 80000000E-03  a 

CARD 

10 

KDUMPN 

m 

2 

KDUMPT 

m 

2 

IFILE1 

m 

IFILE2 

m 

IFILE3 

m 

IFILEI 

m 

CARD 

11 

NOUTF 

m 

1 

NOUTP 

m 

1 

CARD 

12 

IDRAC 

m 

1 

CARD 

13 

TOLSG 

m 

0. 10000000E-02 

TOLEG 

m 

0. 10000000E-02 

1 

TOLUG 

m 

0. 10000000E-02 

j 

TOLPG 

m 

0. 10000000E-02 

TOLRG 

m 

0. 10000000E-02 

j 

CARD 

14 

TOLSP 

m 

0. 10000000E-02 

TOLEP 

m 

Not  used  (an  arbitrary  value  of 

0. 10000000E-02 

is  read  in) 

TOLUP 

m 

0. 10000000E-02 
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